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PREFACE 


Bounded  nominal  paths  can  be  constructed  in  the  vicinity  of  the 
interior  equilibrium  point  (sometimes  called  a  libration  or  Lagrange 
point)  for  the  Sun-Earth+Moon  Elliptic  Restricted  Three-Body  Problem. 
Numerical  integration  is  used  to  generate  the  periodic  or  quasi- 
periodic  reference  trajectories  in  this  effort,  and  the  numerical  data 
is  then  curve  fit  using  a  cubic  spline  routine.  The  force  model  used 
in  this  effort  includes  solar  radiation  pressure,  the  gravitational 
attractions  of  the  Sun  and  the  Earth+Moon  barycenter,  and  the 
centrifugal  force  associated  with  rotation  of  the  system.  A  spacecraft 
near  a  libration  point  orbit  between  the  Earth  and  the  Sun  can  study 
the  interaction  of  the  Sun’s  corona  with  the  terrestrial  environment 
and  will  thus  be  of  great  scientific  value. 

The  spacecraft  will,  however,  drift  from  the  nominal  path,  and  the 
forces  affecting  the  spacecraft  orbit  have  differing  levels  of 
uncertainty.  Both  range  and  range-rate  tracking  also  include 
inaccuracy  in  the  measurement.  The  accumulated  error  in  the 
spacecraft’s  position  and  velocity  relative  to  the  nominal  path  after  a 
predetermined  period  of  tracking  can  be  computed.  This  error,  or 
uncertainty,  in  the  spacecraft  state  is  measured  through  simulations, 
commonly  referred  to  as  orbit  determination  error  analysis,  and  is 
presented  as  either  variances  or  standard  deviations  of  the  state 
vector  elements. 

The  state  uncertainty  computed  in  the  error  analysis  can  then  be 
input  to  a  station-keeping  algorithm.  The  algorithm  computes  control 
manuevers  that  return  the  spacecraft  to  the  vicinity  of  the  nominal 
(unstable)  path.  A  control  algorithm  is  required  for  an  interior 
libration  point  orbit,  and  variations  in  orbital  shapes  and  sizes  may 
have  some  effect  on  the  station-keeping  costs.  Several  algorithms  are 
derived  and  are  used  to  test  differences  in  station-keeping  costs. 

This  effort  is  supported  by  the  Frank  J.  Seiler  Research 
Laboratory  and  has  been  conducted  as  doctoral  research  under  the 
direction  of  Professor  K.C.  Howell,  School  of  Aeronautics  and 
Astronautics,  Purdue  University,  West  Lafayette,  Indiana. 
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INTRODUCTION 


With  the  expansion  of  space  exploration  programs  worldwide, 
interest  has  increased  in  the  design  of  innovative,  complex,  and  yet 
low-cost  spacecraft  trajectories  that  meet  demanding  mission 
requirements.  In  most  of  the  missions  flown  in  the  last  few  decades, 
the  spacecraft  spent  the  majority  of  the  fligM  time  in  a  force 
environment  dominated  by  a  single  gravitational  field.  For  the 
preliminary  mission  analysis  in  these  cases,  additional  attracting 
bodies  and  other  forces  could  be  modeled,  when  required,  as  perturbing 
influences.  Analysis  of  some  recently  proposed  and  more  adventurous 
missions,  such  as  those  involving  libration  point  orbits,  will  require 
dynamic  models  of  higher  complexity,  since  at  least  two  gravitational 
fields  are  of  nearly  equal  influence  on  the  spacecraft  throughout  the 
majority  of  the  mission.  Thus,  trajectories  determined  for  a  system 
consisting  of  numerous  gravitational  forces  have  been  of  particular 
theoretical  and  practical  interest  in  recent  years. 

One  type  of  many-body  problem,  motion  within  a  three-body  system 
of  forces,  has  a  wide  range  of  applications.  The  general  problem  of 
three  bodies  assumes  that  each  body  has  finite  mass  and  that  the  motion 
is  a  result  of  mutual  gravitational  attraction.  When  the  mass  of  one 
of  the  three  bodies  is  assumed  to  be  sufficiently  small  (infinitesimal) 
so  that  it  does  not  affect  the  motion  of  the  other  two  bodies 
(primaries)  in  the  system,  the  "restricted  three-body  problem"  results. 
The  primaries  can  be  further  assumed  to  be  moving  in  known  elliptic  or 
circular  orbits  about  their  common  center  of  mass.  Therefore,  the 
elliptic  restricted  three-body  problem,  where  the  primaries  are  assumed 
to  be  in  known  elliptic  orbits,  may  be  considered  a  reasonably 
approximate  model  for  a  spacecraft  moving  within  the  gravitational 
fields  of  the  Sun  and  the  Earth,  for  instance. 
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In  the  .ormulation  of  the  restricted  three-body  problem,  one  mass 
is  defied  as  infinitesimal  relative  to  the  remaining  two  masses 
(primaries).  The  primaries,  unaffected  by  the  infinitesimal  mass,  move 
under  their  mutual  gravitational  attractions.  In  the  elliptic 
restricted  three-body  problem  (ER3BP),  the  primaries  are  assurr.  1  to 
move  on  elliptic  paths.  If  the  eccentricity  of  the  primaries’  orbit  is 
assumed  to  be  zero,  the  circular  restricted  three-body  problem  (CR3BP) 
results.  Even  for  known  primary  motion,  however,  a  general, 
closed-form  solution  for  motion  of  the  third  body  of  infinitesimal  n.oss 
does  not  exist.  In  the  restricted  three-body  problem  (ER3BP  or  CR3BP), 
five  equilibrium  (libration)  solutions  can  be  found.  These  equi librium 
points,  sometimes  called  Lagrange  points,  are  particular  solutions  of 
the  equations  of  motion  governing  the  path  of  the  infinitesimal  mass 
moving  within  the  gravitational  fields  of  the  primaries. 

The  equilibrium  points  are  defined  relative  to  a  coordinate  system 
rotating  with  the  primaries.  At  these  locations,  the  forces  on  the 
spacecraft  are  in  equilibrium.  These  forces  include  the  gravitational 
forces  from  the  massive  bodies  and  the  centrifugal  force  associated 
with  the  rotation  of  the  system.  (The  addition  of  solar  radiation 
pressure  to  the  force  model  changes  the  locations  of  the  five  Lagrange 
points,  although  they  can  still  be  defined,  and  these  solar  radiation 
effects  are  discussed  in  Gordon111.)  The  libration  points  are  located 
in  the  plane  of  primary  rotation.  Three  of  the  libration  points  are  on 
the  line  between  the  two  massive  bodies,  and  one  of  these  col  linear 
points  is  interior  to  the  primaries.  The  last  two  points  are  at  the 
vertices  of  two  equilateral  triangles  in  the  plane  of  primary  rotation. 
The  triangles  have  a  common  base  that  is  the  line  between  the  primary 
masses . 

For  the  CR3BP,  the  five  libration  points  are  stationary  relative 
to  the  rotating  reference  frame.  If  the  problem  is  generalized  to  the 
ER3BP,  the  libration  points  pulsate  as  the  distance  between  the 
primaries  varies  with  time.  In  both  the  circular  and  elliptic 
restricted  problems,  two-dimensional  and  three-dimensional 
trajectories,  both  periodic  and  quasi-per iodic  paths,  can  be  computed 
in  the  vicinity  of  these  libration  points. 
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Three-dimensional,  periodic  “halo"  orbits  in  the  vicinity  of  the 

collinear  libration  points  have  been  studied  since  the  late  1960s. 

Early  work  concerning  these  orbits  was  motivated  by  studies  related  to 

exploring  the  far  side  of  the  Moon.  These  studies  were  completed  in 

support  of  the  planned  Apollo  18  lunar  exploration  mission  that  was 

later  canceled.  Robert  Farquhar  coined  the  term  "halo”  to  describe  a 

three-dimensional,  periodic  orbit  near  a  libration  point  on  the  far 

[2] 

side  of  the  Moon  in  the  Earth-Moon  system.  These  orbits  were 

designed  to  be  large  enough  so  that  the  spacecraft  would  be  constantly 

in  view  of  the  Earth  and  would  thus  appear  as  a  halo  around  the  Moon. 

A  communications  station  in  this  type  of  orbit  could  maintain  constant 

contact  between  the  Earth  and  a  lunar  experimentation  station  on  the 

f31 

far  side  of  the  Moon. 

Quasi-periodic  orbits  near  libration  points  are  also  currently  of 

great  research  interest.  The  variations  in  size  and  shape  that  a 

quasi-periodic  orbit  can  exhibit  may  add  valuable  flexibility  for 

mission  planning.  This  type  of  bounded,  three-dimensional  libration 

point  trajectory  is  called  a  Lissajous  orbit  since  specific  planar 

projections  of  these  quasi-periodic  trajectories  may  lock  like  a 

special  type  of  "Lissajous"  curve.  Physicist  Jules  Antoine  Lissajous 

(1822-1880)  investigated  curves  that  were  generated  by  compounding 

simple  harmonic  motions  at  right  angles,  and  he  delivered  a  paper  on 

this  subject  to  the  Paris  Academy  of  Sciences  in  1857.  Nathaniel 

Bowditch  of  Salem,  Massachusetts,  had  conducted  some  similar  work  in 

1815.  Lissajous  curves  have  a  wide  variety  of  shapes  that  depend  on 

the  frequency,  phase,  and  amplitude  of  the  orthogonal  components  of  the 
[4  5] 

motion.  ’  When  the  in-plane  and  the  (orthogonal)  out-of-plane 

frequencies  of  the  linearized  motion  are  nearly  (but  not)  equal,  the 

resulting  path  is  typically  called  a  Lissajous  trajectory. 

A  method  to  generate  approximations  for  this  type  of 

quasi-periodic  orbital  path  was  developed  analytically  by  Farquhar  and 

Kamel  in  1973. 161  They  derived  a  third-order  approximate  analytic 

solution  for  a  transiunar  libration  point  orbit  in  the  Earth-Moon  ER3BP 

that  also  included  solar  gravity  perturbations.  In  1975,  Richardson 

and  Cary  then  developed  a  fourth-order  analytic  Lissajous  approximation 

17] 

in  the  Sun-Earth+Moon  barycenter  system.  The  notation  "Earth+Moon" 
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indicates  that  the  Earth  and  the  Moon  are  treated  as  one  body  with  mass 

center  at  the  Earth-Moon  barycenter.  In  consideration  of  the  lunar 

perturbation,  Farquhar  has  shown  that  the  accuracy  of  solutions  in  the 

Sun-Earth  CR3BP  can  be  enhanced  if  the  col linear  libration  points  are 

defined  along  the  line  between  the  Sun  and  the  center  of  mass  o:  *he 

[8] 

Earth  and  the  Moon.  Since  1975,  a  few  researchers  have  considered 

methods  to  numerically  generate  Lissajous  trajectories,  but  the  lack  of 

periodicity  of  a  Lissajous  path  complicates  numerical  construction  of 

bounded  trajectories.  Howell  and  Pernicka  have  developed  a  numerical 

technique  for  determination  of  thrr 3-dimensional ,  bounded  Lissajous 

[9-14] 

trajectories  of  arbitrary  size  and  duration.  Orbits  computed 

with  their  method  are  used  in  this  effort  to  define  the  nominal  path 
near  which  the  spacecraft  will  be  maintained. 

Trajectory  determination  for  a  spacecraft  that  moves  under  the 
influence  of  a  two-body  system  of  forces  will,  however,  be  affected  by 
many  sources  of  error,  including  tracking  errors,  modeling  uncertainty, 
and,  possibly,  control  input  errors.  Orbit  determination  error 
analysis  seeks  to  quantify  the  impact  of  the  numerous  errors  that 

affect  the  motion  of  the  spacecraft.  The  result  of  the  error  analysis 
is  a  determination  of  the  spacecraft  position  and  velocity  uncertainty 
after  some  predetermined  period  of  flight  during  which  the  spacecraft 
is  affected  by  both  the  uncertainties  in  the  forces  and  the  errors  in 
tracking  data.  The  combined  magnitude  of  the  errors  may  be  found  to 
vary  depending  on  the  size  and  shape  of  the  spacecraft  orbit.  A 

reduction  in,  or  a  more  accurate  estimation  of,  the  magnitudes  of  the 

individual  errors  may  be  possible  and  could  then  lead  to  a  significant 

reduction  in  overall  vehicle  position  and  velocity  uncertainty. 

This  reduced  level  of  position  and  velocity  uncertainty  may,  in 
turn,  reduce  orbital  "maintenance"  costs,  such  as  the  propellant 
required  to  keep  the  spacecraft  near  the  nominal  orbit.  The  orbital 
maintenance  routine  is  referred  to  here  as  "station-keeping."  This 
cost  is,  in  part,  dependent  on  the  accuracy  of  the  tracking  information 
because  position  updates  using  inaccurate  tracking  data  may  result  in 
inefficient  use  of  control  energy  and  may  also  lead  to  unacceptable 
spacecraft  drift  from  the  nominal  path.  Other  error  sources  may  also 
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affect  spacecraft  drift  from  the  (unstable)  reference  trajectory  and, 
therefore,  may  increase  station-keeping  costs. 

This  research  is  concerned  with  developing  and  evaluating 
station-keeping  algorithms  for  libration  point  orbits.  The  errors 
derived  in  orbit  determination  error  analysis  studies  are  used  as 
random  inputs  in  Monte  Carlo  simulations  of  the  competing 
station-keeping  algorithms.  Other  random  inputs  include  solar 
radiation  pressure  uncertainty  and  control  input  errors.  The  output  of 
the  station-keeping  trials  is  a  function  of  the  several  random  inputs 
and  is  consequently  treated  as  a  random  variable.  Statistical  goodness 
of  fit  tests  and  equivalence  of  means  and  variances  tests  can  then  be 
appropriately  conducted.  The  results  associated  with  orbits  designed 
to  be  periodic  and  quasi-periodic  orbits  are  also  compared.  Chapter  1 
briefly  summarizes  the  background  of  the  elliptic  restricted  three-body 
problem.  Chapter  2  then  derives  several  station-keeping  methods,  and, 
finally,  Chapter  3  covers  the  results  obtained. 
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CHAPTER  1:  BACKGROUND 


In  this  chapter,  the  elliptic  restricted  three-body  problem  and 
the  associated  coordinate  systems  are  reviewed;  the  equations  of  motion 
for  an  infinitesimal  mass  moving  in  the  gravity  fields  of  two  massive 
bodies  are  then  presented.  Next,  locations  of  the  libration  points  are 
discussed.  The  state  transition  matrix  and  the  construction  of  bounded 
nominal  orbits  near  the  col linear  Lagrange  points  are  then  summarized. 
Finally,  curve  fitting  the  nominal  trajectory  is  covered.  A  more 
thorough  discussion  of  these  topics  is  presented  in  Gordon. 11,155 


A.  Elliptic  Restricted  Three-Body  Problem 

The  elliptic  restricted  three-body  problem  is  a  simplification  of 
the  general  problem  of  three  bodies.  In  the  general  three-body 
problem,  each  of  the  three  bodies  is  assumed  to  be  a  particle  of  finite 
mass  and,  thus,  exerts  an  Influence  on  the  motion  of  each  of  the  other 
bodies.  Neither  the  general  nor  the  restricted  problem  of  three  bodies 
has  a  general  closed-form  solution.  However,  when  problem 
simplifications  are  made,  particular  solutions  can  be  determined.  If 
the  mass  of  one  of  the  bodies  is  restricted  to  be  infinitesimal,  such 
that  it  does  not  affect  the  motion  of  the  other  two  massive  bodies 
(primaries),  the  restricted  three-body  model  results.  The  primaries 
are  assumed  to  be  in  known  elliptic  (ER3BP)  or  circular  (CR3BP)  orbits 
about  their  common  mass  center  (barycenter ) .  The  problem  can  then  be 
completely  described  by  a  single  second-order  vector  differential 
equation  with  variables  appropriately  defined  for  a  specified 
coordinate  frame. 
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B.  Coordinate  Systems 


The  two  standard  coordinate  systems  used  in  the  analysis  of  this 

problem  have  a  common  origin  at  the  center  of  mass  (barycenter)  of  the 

primaries.  Primaries  with  masses  m  and  m  such  that  m  a  m  are 

12  12 

assumed  here,  although  this  distinction  is  arbitrary.  The 

infinitesimal  mass  is  denoted  as  m  .  These  masses  (m  ,m  ) 

3  12  3 

correspond  to  particles  situated  at  points  P  ,  P  ,  and  P  , 
respectively.  The  barycenter  is  denoted  as  "B,"  and  the  resulting 
arrangement  is  shown  in  Figure  1-1.  The  rotating  coordinate  system  is 
defined  as  x  y  z  ,  and  the  Inertial  system  is  identified  as  x  y  z  . 

R  R  R  III 

Note  that  both  coordinate  systems  are  right-handed,  and  the  x  and  y 

axes  for  both  systems  are  in  the  plane  of  motion  of  the  primaries.  The 

xt  axis  is,  of  course,  assumed  to  be  oriented  in  some  fixed  direction; 

in  this  specific  formulation  of  the  problem,  it  is  assumed  to  be 

parallel  to  a  vector  defined  with  a  base  point  at  the  Sun  and  directed 

toward  periapsis  of  the  Earth’s  orbit.  The  rotating  x  axis  is  defined 

R 

along  the  line  that  joins  the  primaries  and  is  directed  from  the  larger 
toward  the  smaller  primary.  The  z  axes  are  coincident  and  are  directed 
parallel  to  the  primary  system  angular  momentum  vector.  The  y  and  y 

I  R 

axes  complete  the  right-handed  x  y  z  and  x  y  z  systems,  respectively. 

Ill  R  R  R 
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Pigure  1-1.  Coordinate  Systems  With  Barycenter  Origin. 
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C.  Equations  of  Notion 


Newtonian  mechanics  are  used  to  formulate  the  equations  of  motion 
for  m3  (the  spacecraft)  relative  to  B  as  observed  in  the  inertial 
reference  frame.  The  sum  of  the  forces  on  m3  resulting  from  both  the 
gravity  fields  of  masses  (the  Sun)  and  m2  (the  Earth-Moon 
barycenter)  and  from  the  solar  radiation  pressure  can  be  used  to 
produce  the  following  second-order  vector  differential  equation: 


It 

P 


m 


=  -  G  (- 


l  *  ? 


m 


)  d  -  G  (- 


2  v  - 


)  r  +  ( 


*§-)  d. 


(1-1) 


The  overbar  denotes  a  vector,  and  primes  indicate  differentiation 
with  respect  to  dimensional  time.  All  quantities  are  dimensional,  as 
appropriate,  and  the  quantity  "G"  is  the  universal  gravitiational 
constant.  The  scalars  "d"  and  "r"  in  equation  (1-1)  denote  the 
magnitudes  of  the  vectors  d  and  r,  respectively,  as  depicted  in 
Figure  1-1.  The  dimensionless  scalar  "k"  is  the  solar  reflectivity 
constant,  and  "S"  is  the  solar  radiation  pressure  constant.  The 
formulation  of  the  solar  radiation  force  model  and  the  values  for  the 
solar  radiation  constants  are  derived  from  previous  work  by  Bell.1161 
The  values  of  the  constants  are  described  in  Gordon. [l1 

The  position  vector  p  is  written  in  rotating  components  as 


p  =  x  x  +  y  y  +  zz 

R  JR  R 


(1-2) 


where  x  ,y  ,z  are  unit  vectors.  The  velocity  and  the  acceleration  of 

R  R  R 

the  spacecraft  (particle  P3  with  mass  m3)  relative  to  the  barycenter  B 
as  observed  in  the  inertial  reference  frame  can  then  be  described.  The 
following  kinematic  expression  for  p"  can  be  derived: 
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(1-3) 


p"  =  (x"-0"y-20'y,-0'2x)x  +(y"+0"x+20'x' -0,2y)y  +z"z  . 

R  R  R 


Three  scaled  equations  of  motion  for  P3  can  be  derived  using  the 
the  following  scaling  factors: 


(1)  The  sum  of  the  masses  of  the  primaries  equals  one 
mass  unit. 

(2)  The  mean  distance  between  the  primaries  equals  one 
unit  of  distance. 

(3)  The  universal  gravitational  constant  is  equal  to  one 
unit  by  proper  choice  of  characteristic  time. 


The  dimensional  equations  of  motion  can  be  simplified  and  scaled 
by  introducing  the  characteristic  quantities  defined  above  and  by 
introducing  the  nondimensional  mass  ratio  p,  "psuedo-potential"  U,  and 
the  scaled  solar  radiation  constant  s: 


and 


m  +  m 
l  z 


(1-4) 


U  = 


(1-p) 

d 


+ 


£ 

r 


•  p  2  2 

0  (x  +  y 


k  s 
d 


(1-5) 


where  the  dot  denotes  the  derivative  with  respect  to  characteristic 
time.  The  scaled  solar  radiation  constant,  s,  is  derived  by  using  the 
characteristic  quantities  described  above.  Then,  the  vector  magnitudes, 
"d"  and  "r,"  are  written  in  terms  of  scaled  quantities  as: 
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(1-6) 


,  ii  „,2  2  2,1/2 

d  =  [  (x  +  |i  R)  +  y  +  z  ]  , 


II  r,  n,2  2  2,1/2 

r  =  ( (x  -  R  +  (i  R)  +  y  +  z] 


(1-7) 


The  three  scalar  second-order  differential  equations  that  result 
can  be  written  in  terms  of  characteristic  quantities  as 


x  -  2  9  y  = 


0  y  =  U  +  0  y, 

X 


(1-8) 


y  +  2  0  x  = 


-  0  x 


■  U  -  0  x, 
y 


(1-9) 


3U 

di 


(1-10) 


If  the  primaries  are  assumed  to  be  moving  in  a  circular  orbit, 
equations  (1-8),  (1-9),  and  (1-10)  reduce  to  three  scalar  equations  in 
the  simplified  form: 
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x  -  2  y  = 


au 

5 x 


=  u 


(1-11) 


*•  • 
y  +  2  x 


(1-12) 


(1-13) 


The  scalar  equations  can  be  used  to  locate  the  five  libration 
points  in  the  rotating  reference  frame. 


D.  Locations  of  the  Lagrangian  Points 

By  using  scalar  equations  (1-11),  (1-12),  and  (1-13)  for  motion  in 
the  CR3BP,  the  locations  of  the  stationary  equilibrium  points  can  be 
determined.  Equations  (1-8),  (1-9),  and  (1-10)  can  be  used  to 
determine  ratios  of  distances  that  are  constant  in  the  ER3BP;  these 
ratios  are  related  to  the  locations  of  libration  points  that  have  been 
defined  in  the  ER3BP  and  that  "pulsate"  with  respect  to  the  rotating 
reference  frame  as  the  distance  between  the  primaries  varies  with  time. 


1 .  The  CR3BP 

In  the  CR3BP,  the  five  libration  points  are  equilibrium  points  and 
are  stationary  with  respect  to  the  rotating  coordinate  frame,  that  is, 
they  are  locations  at  which  the  forces  on  the  third  body  sum  to  zero. 
The  arrangement  of  points  and  the  corresponding  nondimensional 
distances  are  depicted  in  Figure  1-2.  Note  that  three  of  the  libration 
points  (L^  L2>  L3)  are  collinear  with  the  primaries;  one  collinear 


12 


point  (L  )  is  interior  to  the  primaries.  The  remaining  two  points  (L 
and  Lg)  are  located  at  the  vertices  of  two  equilateral  triangles  that 
are  in  the  plane  of  primary  rotation  and  that  have  a  common  base 
between  the  primaries. 

In  the  CR3BP,  the  libration  points  are  stationary  in  the  rotating 
coordinate  frame.  Stationary  points  are  defined  as  points  at  which  the 
relative  velocity  and  acceleration  are  zero,  such  that 

x  =  y  =  z  =  x  =  y  =  z  =  0.  (1-14) 


By  using  equations  (1-14)  in  equations  (1-11)  through  (1-13),  the 
useful  conditions  U  =  U  =  U  =  0  are  found.  The  three  collinear 

x  y  z 

libration  points  can  be  readily  located  by  further  noting  that 
y  =  z  =0  for  the  points  located  on  the  rotating  xr  axis. 


2.  The  ER3BP 

Five  libration  points  also  exist  in  the  ER3BP,  but  they  are  not 
stationary  relative  to  the  rotating  frame;  rather,  the  collinear  points 
pulsate  along  the  x  axis,  and  the  triangular  points  pulsate  relative 

K 

to  both  the  x  and  the  y  axes  as  the  distance  between  the  primaries 

R  R 

varies  with  time.  The  equilibrium  solutions  can  be  located  by  using 
equations  (1-8)  through  (1-10)  to  find  ratios  of  certain  distances 
that  are,  in  fact,  constant  in  the  problem.  The  collinear  libration 
points  in  the  ER3BP  can  be  found  by  assuming  x  *  0,  x  *  0,  and  y  =  y  = 
z  =  z  =  y  =  z  =  0.  The  relative  locations  of  the  libration  points  in 
the  ER3BP  are  depicted  in  Figure  1-3. 
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Figure  1-3.  Lagrange  Point  Locations  in  the  Scaled  ER3BP. 


E.  State  Transition  Matrix 


The  state  transition  matrix  is  used  in  the  calculation  of  the 
acceptable  nominal  trajectory,  and  it  must  also  be  available  at  varying 
time  intervals  along  the  nominal  path  for  orbit  determination  error 
analysis  investigations  and  station-keeping  studies.  The  transition 
matrix  is  derived  in  connection  with  a  linearizing  analysis. 

The  equations  of  motion  for  the  Infinitesimal  mass  in  the  ER3BP 
can  be  linearized  about  a  reference  trajectory  (nominal  path)  that  is  a 
solution  of  the  differential  equations.  The  states,  three  position  and 
three  velocity,  and  the  state  vector  x  are  defined  as 


x  =  x,  x  =  y,  x  =z,x  =  x,  x  =  y ,  x  =  z,  (1-15) 


and 


=  lx. 


X.) 


(1-16) 


The  reference  trajectory  is  defined  as  x  .  Therefore,  using  a 
Taylor’s  series  approach,  the  expansion  about  the  reference  path  is 
written  in  the  form  of  the  first-order  variational  equation 


—  (  x)  =  x  =  A( t )  x  (1-17) 

dt 


where  x  =  x  -  x  is  understood  to  be  the  vector  of  residuals  relative 
REF 

to  the  nominal  solution,  and  the  matrix  A ( t )  contains  the  first-order 
terms  in  the  Taylor’s  series  expansion  of  the  equations  of  motion  about 
the  nominal  or  reference  solution  of  interest. 
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Using  equations  (1-8)  through  (1-10),  Aft)  can  be  expressed  as 


Aft) 


0  I 

Urr+0(2  29(2 


(1-18) 


where  all  four  submatrices  are  dimension  3x3  and 


Uxx 

Uxy 

Uxz 

Uyx 

Uyy 

Uyz 

Uzx 

Uzy 

Uzz 

with 


£2  = 


'  0 
-1 
0 


1 

0 

0 


0 

0 

0 


In  equation  (1-19),  the  notation  is  simplified  for  the  partial 
derivatives;  for  instance 


32U 


dx 


2 


Uxx . 


The  matrix  A(t)  can  then  be  evaluated  along  the  reference  trajectory. 

The  vector  differential  equation  (1-17)  governing  the  state 
variations  from  the  nominal  path  has  a  solution  of  the  form 
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(1-20) 


x ( t )  =  #(t, t  )  x(t  ) 
0  0 


where  $(t,t  )  Is  the  state  transition  matrix  at  time  "t"  relative  to 
0 

time  "  t  . "  The  matrix  ♦  ,  then,  represents  the  sensitivities  of  the 
states  at  time  "t"  to  small  changes  in  the  initial  conditions.  It  is 
determined  by  numerically  integrating  the  matrix  differential  equation 


-r  ♦ft-1  )  =  *(t,t  )  =  A ( t )  *(t,t  ), 

dt  0  0  0 


(1-21) 


with  initial  conditions  $(t  ,  t  )  =  I,  the  6x6  identity  matrix.  Thus, 

0  0 

the  nonlinear  equations  of  motion  in  (1-8)  through  (1-10)  and  the 
matrix  equation  (1-21)  combine  to  result  in  42  first-order  differential 
equations  that  can  be  simultaneously  integrated  numerically  to 
determine  the  state  vector  and  its  associated  transition  matrix  at  any 
instant  of  time.  The  reference  trajectories  that  are  of  interest  in 
this  research  are  generated  by  a  numerical  integration  method  that  uses 
a  differential  corrections  process  developed  by  Howell  and 
Pernicka.  19  141  The  orbits  include  solar  radiation  pressure  forces  as 
formulated  by  Bell1161  specifically  for  an  orbit  associated  with  the 
interior  Lagrange  point  in  the  Sun-Earth  system.  The  numerical 
integration  routines  used  in  this  work  are  fourth-  and  fifth-order 

f  17  j 

Runge-Kutta  formulas  available  in  the  386-Matlab  software  package. 


F.  Bounded  Orbits  Near  Libration  Points 

The  computation  of  bounded  periodic  and  quasi-periodic  orbits  in 
the  vicinity  of  libration  points  has  been  of  increasing  interest  during 
the  past  100  years.  This  section  first  discusses  the  stability  of  the 
libration  points  in  the  CR3BP  and  the  ER3BP.  The  construction  of 
bounded  orbits  near  the  collinear  Lagrange  points  is  then  summarized. 
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Finally,  the  specific  reference  trajectories  used  In  the  orbit 
determination  error  analysis  and  station-keeping  studies  in  this  work 
are  introduced. 


1.  Stability  of  the  Libration  Points  in  the  CR3BP 

The  accomplishments  of  those  researchers  who  have  constructed 

bounded  orbits  near  collinear  libration  points  are  particularly 

significant  because  the  collinear  points  are  considered  "unstable" 

points  of  equilibrium  but  with  (only)  one  mode  producing  positive 

exponential  growth.  Bounded  motion  in  their  vicinity,  therefore,  is 

determined  by  deliberately  not  exciting  the  unstable  mode.  A  second 

mode  produces  negative  exponential  orbital  decay  and  is  also 

deliberately  not  excited.  In  the  CR3BP,  the  remaining  four  eigenvalues 

are  purely  imaginary.  The  existence  of  initial  conditions  that  result 

in  only  trigonometric  (sinusoidal)  functions  as  solutions  means  that 

the  collinear  libration  points,  while  unstable,  possess  conditional 

stability  (with  proper  choice  of  initial  conditions)  in  the  linear 
[18] 

sense. 

The  triangular  libration  points  are  marginally  stable  in  the 

linear  sense  for  a  specific  range  of  primary  mass  ratio  in  the  CR3BP. 

Purely  imaginary  roots  in  two  conjugate  pairs  are  obtained  for  jis.0385, 

which  is  given  here  to  four  decimal  places  and  is  sometimes  referred  to 
( 19  ] 

as  Routh’ s  value.  The  mass  ratios  (listed  here  to  three  decimal 

places),  for  example,  in  the  three-body  systems  of  the  Earth-Moon 
(p  =  1.216  x  10  6),  Sun-Ear th+Moon  (ji  =  3.022  x  10  6)  and  Sun-Jupiter 

-4 

(p  =  9.485  x  10  )  all  satisfy  the  mass  ratio  requirement  for  marginal 

stability  of  the  triangular  points  in  the  linearized  model.  Natural 
satellites,  such  as  the  Trojan  asteroids  or  a  moon  of  Saturn,  occupy 
linearly  stable  orbits  near  triangular  libration  points  in  their 
respective  systems. 
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2.  Stability  of  the  Libration  Points  in  the  ER3BP 


Several  researchers  have  analyzed  the  stability  of  the  libration 

points  in  the  elliptic  problem,  where  both  the  mass  ratio,  p,  and  the 

[18-22] 

primary  orbit  eccentricity,  e,  influence  stability.  The 

instability  of  the  collinear  libration  points  as  determined  in  the 

circular  problem  for  all  the  values  of  mass  parameter  persists  for  the 

elliptic  problem;  an  analysis  of  the  collinear  points  shows  instability 

for  any  combination  of  the  values  of  both  p  and  e. 

The  results  of  a  linearized  stability  analysis  regarding  the 

effects  of  eccentricity  and  mass  ratio  on  the  linear  stability  of  the 

121] 

triangular  points  have  been  published  by  Danby  and  then  later  by 
[22] 

Bennett  .  Both  Danby  and  Bennett  have  numerically  generated  graphic 
depictions  of  the  linear  stability  region  in  the  p-e  plane.  For  the 
eccentricity  in  the  Sun-Earth+Moon  ER3BP,  the  value  of  p  which  ensures 
linear  stability  is  only  slightly  less  than  Routh’s  value  (decreased  by 
approximately  one  percent).  An  interesting  aspect  of  the  p-e  stability 
region  is  thaL  a  range  of  values  of  p  greater  than  Routh’s  value  also 
defines  a  region  of  linear  stability  for  a  specific  range  of  values  of 
e  less  than  .3143. 

3.  Construction  of  Bounded  Collinear  Libration  Point  Orbits 

The  initial  goal  in  the  process  of  generating  bounded  orbits  near 
a  collinear  (unstable)  libration  point  is  to  avoid  exciting  the 
unstable  mode  associated  with  the  linearized  motion.  The  meteoric  dust 
particles  that  may  be  orbiting  near  Lagrange  point  L in  the  Sun-Earth 
system  could  only  linger  near  that  point  if  they  arrive  with  the 
"correct"  initial  position  and  velocity  states  relative  to  L  .  The 
"correct"  initial  conditions  will  only  (primarily)  excite  the  stable 
modes  associated  with  the  linearized  motion  and  not  (or  minimally) 
excite  the  unstable  mode.  The  degree  to  which  the  unstable  mode  is 
excited  will  determine  the  length  of  time  that  the  dust  particles 
linger  near  L . 
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The  third-order  analytic  representation  is  used  in  this  work  to 
provide  the  initial  model  for  the  trajectories.  The  method  of 
successive  approximations,  using  the  linearized  solution  as  the  first 
approximation  to  the  nonlinear  orbital  path,  and  the  method  of  dual 

1 6  7  23 ] 

time  scales  are  used  to  derive  the  third-order  result.  ’  ’  The 

method  of  successive  approximations  is  used  to  generate  an  asymptotic 
series  in  an  appropriately  small  parameter.  (The  square  root  of  the 
eccentricity  of  the  primary  orbit,  that  is  the  orbit  of  the  Earth-Moon 
barycenter  about  the  Sun,  is  the  small  parameter  used  here. )  The 
method  of  dual  time  scales  is  used  to  convert  the  system  of  ordinary 
differential  equations  into  a  system  of  partial  differential  equations. 
In  general,  the  method  of  multiple  scales  allows  the  various  nonlinear 
resonance  phenomena  to  be  included  in  the  approximate  analytic  solution 
and  provides  a  method  to  remove  secular  terms.  (Here,  "secular"  refers 
to  terms  that  include  the  time  variable  and  is  derived  from  the  French 
"siecle"  meaning  century. ) 

[7J 

The  analytic  solution  of  Richardson  and  Cary  for  the 
Sun-Ear th+Moon  ER3BP  has  been  derived  to  fourth  order,  but  the  third- 

19-14] 

order  approximation  is  found  to  be  sufficient  for  this  research. 

A  numerical  Integration  algorithm,  using  a  differential  corrections 
procedure  that  is  designed  to  adjust  the  first  guess  as  obtained  from 
the  analytic  approximation,  can  then  be  used  to  numerically  generate 

[9-14] 

the  orbit  of  interest.  A  method  developed  by  Howell  and  Pernicka 
is  used  here  to  generate  the  orbital  paths.  Their  method  initially 
employs  the  approximate  analytic  solution  to  compute  target  points.  A 
two-level  (position  matching  then  velocity  matching),  multi-step 
differential  corrections  algorithm  is  used  to  construct  a  numerically 
Integrated,  bounded  trajectory  that  is  continuous  in  position  and 
velocity.  A  solar  radiation  pressure  model  developed  by  Bell1161  is 
also  incorporated  in  the  numerical  integration  procedure. 

The  method  of  Howell  and  Pernicka,  including  solar  radiation 
pressure,  uses  an  initial  analytic  guess  that  represents  a  halo  orbit 
or,  alternatively,  a  considerably  smaller  Lissajous  path.  The 
higher-order  terms  tend  to  slightly  alter  the  first-order  periodic  or 
quasi-periodic  path.  Consequently,  the  initial  target  path  for  a  halo 
orbit  will  generally  not  be  precisely  periodic.  The  two-level, 
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multi-step  differential  corrections  procedure  then  adjusts  the  initial 
analytic  target  orbit  and,  therefore,  will  compute  a  halo-type  orbit 
that  is  nearly  (but  not  exactly)  periodic. 


4.  The  Reference  Paths  Generated  for  This  Work 

Precisely  periodic  halo  orbits  exist  in  the  CR3BP.  They  also 
exist  in  the  ER3BP,  but,  in  the  ER3BP,  they  are  multiple  revolution 
trajectories  with  periods  much  longer  than  those  of  interest  here. 
Nearly  periodic  orbits  are  more  practical  in  the  ER3BP  and  are  much 
more  likely  to  be  used  in  mission  planning;  therefore,  the  goal  here 
should  be  slightly  modified  to  be  the  comparison  of  Lissajous  and 
"halo-type"  orbits.  The  general  shapes  of  the  three-dimensional 
halo-type  and  Lissajous  orbits  can  be  seen  by  plotting  three 
orthographic  views  of  each  orbit,  using  the  tabular  data  from  the 
numerical  integration  routine.  Figure  1-4  depicts  three  orthographic 
views  of  point  plots  for  the  Lissajous  orbit  used  in  this  research. 
Figure  1-5  contains  three  orthographic  views  (on  a  slightly  different 
scale)  of  the  considerably  larger  halo-type  orbit.  (Note  that,  in 
general,  the  amplitude  ratio  for  Lissajous  trajectories  is  arbitrary. 
In  halo  orbits,  however,  constraining  the  amplitude  ratio  results  in 
equalized  frequencies  for  in-plane  and  out-of -plane  motion. )  The 
orbits  are  depicted  in  the  rotating  reference  frame  centered  at  L  . 

Both  orbits  are  clearly  not  periodic;  a  Lissajous  orbit  is  often 
called  a  quasi-periodic  path,  and  these  two  orbits  could  clearly  be 
termed  quasi-periodic  or  Lissajous  paths.  The  major  difference  between 
the  orbits  is  the  larger  size  of  the  halo-type  orbit;  however,  other 
differences  are  also  present.  The  maximum  x  and  y  excursions  of  the 
halo-type  orbit  are  approximately  four  times  as  large  as  those  of  the 
Lissajous  path.  Furthermore,  the  direction  of  motion  (clockwise  versus 
counterclockwise),  as  viewed  in  the  y-z  orthographic  depiction,  is 
different  for  the  two  orbits  used  here.  The  direction  of  motion  on  the 
halo-type  orbit  is  counterclockwise  in  the  y-z  depiction;  the  direction 
of  motion  is  clockwise  in  the  y-z  depiction  for  the  Lissajous  path. 
(Both  orbits  include  clockwise  motion  in  the  x-y  depiction. ) 
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indicates  direction  of  motion. 

.  indicates  the  Iteration  point  location. 


Figure  1-5.  Three  Orthographic  Views  of  a  Halo-Type  Orbit. 
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The  two  orbits  can  also  be  differentiated  in  terms  of  the 
direction  of  the  maximum  z  excursion  in  the  x-z  depiction.  If  the 
maximum  z  excursion  is  in  the  positive  z  direction,  the  orbit  can  be 
termed  a  member  of  a  "northern  family"  of  orbits.  When  the  maximum  z 
excursion  of  the  orbit  is  in  the  negative  z  direction,  the  orbit  is 
termed  a  member  of  a  "southern  family"  of  orbits.  In  the  x-z 
orthographic  depiction,  the  smaller  (Lissajous)  path  can  be  seen  to  be 
a  member  of  a  northern  family  of  orbits,  while  the  halo-type  orbit  is  a 
member  of  a  southern  family  of  orbits. 

Future  work  with  these  two  orbits  will  include  studies  that 
generally  require  access  to  a  nominal  path  that  is  at  least  piecewise 
smooth.  Some  method  of  curve  fitting  the  numerically  integrated  data 
must  consequently  be  investigated. 


5.  Curve  Fitting  the  Nominal  Path 

The  primary  goal  of  this  work  is  the  completion  of  orbit 
determination  error  analysis  for  libratlon  point  trajectories.  The 
conventional  method  for  solution  of  state  estimation  problems,  and  the 
technique  used  in  this  effort,  involves  linearizing  the  nonlinear 
equations  of  motion  about  a  reference  solution  (nominal  path)  and  then 
applying  linear  estimation  techniques.  The  orbit  determination  process 
is  thus  changed  from  estimating  the  state  of  a  nonlinear  system  to 
estimating  the  linear,  time-varying  deviations  from  the  reference 
trajectory. 

The  reference  solution  used  in  this  research  is  generated  by 

numerical  integration  of  the  nonlinear  equations  of  motion.  In  one 

study,  an  investigation  that  used  a  consistent  dynamic  model  for  all 

[8] 

comparisons,  Richardson  has  shown  that  a  slight  reduction  in  fuel 
expenditure  can  be  realized  if  a  numerically  integrated,  rather  than  an 
approximate  analytic,  nominal  path  is  computed.  The  numerical 

(9-14] 

integration  method  developed  by  Howell  and  Pernicka  is  used  here 

to  generate  a  set  of  reference  points  for  both  position  (three  states) 
and  velocity  (three  states),  relative  to  the  libration  point  of 
interest,  at  specified  times.  Time  series  point  plots  of  all  six  state 
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variables  for  approximately  a  2-year  segment  of  a  Llssajous  orbit  are 
depicted  in  Figures  1-6  (position  states)  and  1-7  (velocity  states). 
The  method  computes  numerical  data  for  the  six  states  in  a  reference 
frame  that  is  centered  at  the  llbration  point  (in  this  case  L  )  and 
that  rotates  with  the  primaries.  However,  the  state  estimation 
techniques  and  station-keeping  algorithms  used  in  this  work  require 
access  to  a  continuous  nominal  path,  rather  than  point  plots,  of 
acceptable  accuracy. 

The  reference  trajectory,  represented  as  a  (piecewise)  smooth 
curve,  could  be  constructed,  approximately  or  exactly,  through  the 
points  obtained  from  the  numerical  integration  routine.  The  work  here 
assumes  that  a  curve  that  passes  through  the  numerical  aata  (exactly) 
is  preferred.  The  effort  required  to  generate  a  numerical  solution, 
including  forces  modeled  consistent  with  the  ER3BP  (or  even  more 
accurately  modeled  with  ephemeris  data)  would  seem  to  be  wasted  if  the 
reference  curve  deviates  too  far  from  the  numerical  data.  However,  a 
method  that  approximates  a  smooth  curve  through  the  points  is  also 
desirable;  that  is,  linear  interpolation  between  the  numerical  data 

[9] 

points  was  not  considered  acceptable.  In  one  study,  Pernicka  found 
that  station-keeping  costs  for  a  llbration  point  orbit  were,  in  fact, 
sensitive  to  the  accuracy  of  the  curve  fit.  Clearly,  a  piecewise 
linear  curve  fit  could  not  accurately  match  the  concavity  of  the  actual 
orbital  path  between  data  points,  regardless  of  the  size  of  the  time 
steps  used  in  the  numerical  integration  routine. 

Four  methods  of  generating  a  curve  for  the  nominal  trajectory  have 
been  evaluated:  Fourier  series,  least  squares,  weighted  least  squares, 
and  cubic  splines.  The  states  associated  with  a  quasi-periodic  path 
were  thought  to  be  the  most  difficult  to  curve  fit;  therefore,  various 
Lissajous  trajectories  were  used  to  evaluate  the  curve-fitting  methods. 
After  several  curve  fitting  evaluations,  a  cubic  spline  interpolation 
routine  was  selected  to  be  used  to  model  the  reference  trajectory  in 
the  state  estimation  simulations.  The  comparisons  of  the  curve  fit 
methods  are  fully  described  in  Gordon. [1S1 
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Figure  1-6.  Time  Series  Plots  ot  Three  Lissajous  Position  States. 


1000 


27 


G.  Orbit  Determination  Error  Analysis  Results 


Complete,  exact  knowledge  of  the  state  of  a  spacecraft  in  orbit  is 
generally  not  possible.  Individual  state  variables  cannot  be  measured 
precisely,  and  available  measurements  are  usually  some  function  of 
these  state  variables.  For  instance,  a  spacecraft  moving  along  a 
libration  point  trajectory  in  the  Sun-Earth  system  may  be  tracked  using 
range  and  range-rate  measurements  containing  random  errors.  The 
spacecraft  may  be  affected  by  forces  not  included  (or  inadequately 
represented)  in  the  dynamic  model,  and  model  parameters  may  be 
uncertain.  By  definition,  the  linearized  system  of  equations  used  to 
model  the  nonlinear  state  variations  is  a  further  approximation.  Also, 
actual  control  inputs  may  vary  slightly  in  magnitude  and  direction  from 
those  commanded.  These  sources  of  error  make  knowledge  of  the 
spacecraft  state  uncertain.  Computation  of  the  most  likely  current 
state  of  the  spacecraft  in  the  presence  of  measurement  and  model 
uncertainty  is  the  focus  of  orbit  determination. 

Error  analysis  involves  an  investigation  of  the  impact  of  various 
sources  of  error  on  orbit  determination.  The  output  of  an  error 
analysis,  as  used  in  this  work,  provides  the  magnitudes  of  state  vector 
variances  and  covariances,  thus  quantifying  the  relative  contributions 
of  the  significant  error  sources.  This  output  could  then  be  used  to 
predict  how  an  Improvement  in  measurement  accuracy,  for  instance,  would 
lessen  state  uncertainty.  One  benefit  of  more  accurate  knowledge  of 
the  state  might  be  a  reduction  in  station-keeping  costs.  A 
mathematical  procedure  can  be  developed  to  combine  all  information 
about  the  spacecraft  state  and  filter  this  observed  data  based  on  the 
varying  degrees  of  uncertainty.  The  filter  then  produces  a  "best 
estimate"  of  the  state  and  additionally  quantifies  the  resulting  state 
variable  uncertainties  in  preparation  for  an  error  analysis. 

The  orbit  determination  error  analysis  results  can,  in  turn,  be 
used  in  Monte  Carlo  simulations  of  competing  station-keeping 
algorithms.  The  state  error  levels  can  be  modeled  as  random  errors 
with  a  specified  mean  and  probability  distribution.  Of  course,  other 
error  sources,  such  as  solar  radiation  pressure  uncertainty  and  control 


29 


input  errors,  can  also  be  modeled  in  the  Monte  Carlo  simulations. 
[24] 

Gordon  discusses  the  evaluation  of  several  orbit  determination 
error  analysis  methods,  uses  hypothesis  tests  to  determine  the  means 
and  probability  distributions  of  the  errors,  and  computes  the  error 
levels  appropriate  for  use  in  the  station-keeping  simulations.  The 
state  uncertainty  levels  are  found  to  closely  follow  a  zero-mean 
Gaussian  distribution,  and  these  state  error  levels  are  presented  in 
Table  1  in  terms  of  state  element  standard  deviations. 


Table  1.  Error  Levels  Produced  from  Error  Analysis  Studies. 


One  Standard 

Deviation  Levels 

Coordinate 

Halo-Type  Orbit 

Lissajous  Orbit 

x  (km) 

1.46 

1.25 

y  ( km ) 

2.64 

3.35 

z  (km) 

4.81 

3.  19 

x  (mm/sec) 

1.40 

1.25 

y  (mm/sec) 

1.85 

1.41 

z  (mm/sec) 

2.49 

2.51 

CHAPTER  TWO:  STATION-KEEPING  ALGORITHMS 


Because  of  unavoidable  errors,  similar  to  those  associated  with  any 
space  mission,  a  station-keeping  strategy  is  necessary  to  maintain  the 
spacecraft  within  some  predefined  torus  about  the  planned  nominal  path. 
The  size  and  shape  of  the  torus  are  determined  by  mission  objectives 
including  the  possibly  related  requirements  for  orbit  determination, 
scientific  experimentation,  and  minimum  fuel  expenditr  <?.  For  example, 
the  nominal  path  may  be  computed  to  the  highest  degree  of  accuracy  to 
meet  both  scientific  specifications  and  tracking  considerations.  The 
size  of  the  torus  may  then  be  tailored  from  "tight"  to  "loose" 
depending  on  mission  objectives  and  comparison  of  fuel  expenditures 
for  various  options. 

This  chapter  is  organized  into  six  sections.  The  first  section 

contains  an  overview  of  the  general  station-keeping  problem.  The 

second  section  summarizes  the  derivations  of  two  similar  controllers, 

one  of  which  has  been  used  for  libration  point  orbit  station-keeping 

simulations.  The  next  two  sections  summarize  the  derivations  of  two 

related  control  algorithms  formulated  for  this  work.  Both  of  these 

station-keeping  methods  in  the  third  and  fourth  sections  are 

[  2S  2  6] 

innovations  on  methods  previously  derived  by  Dwivedi  ’  and 

[9] 

Pernlcka  The  fifth  section  contains  a  description  of  an  on/off 

controller  developed  for  this  work. 

The  sixth  and  final  section  contains  comparisons  of  the  control 
costs  derived  from  use  of  the  station-keeping  methods  developed  in  this 
work  for  both  a  halo-type  orbit  as  well  as  a  more  general  Lissajous 
orbit.  The  comparisons  require  the  use  of  statistical  hypothesis  tests 
similar  to  those  used  in  the  last  section  of  chapter  two.  A 
station-keeping  simulation  will  produce  a  scalar  value  for  the  total 
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propellant  expended  ( Avt ) ;  also,  each  separate  simulation  will  be 
subject  to  several  random  Inputs  and  can,  therefore,  be  considered  a 
random  trial.  A  group  of  random  trials  (station-keeping  simulations) 
using  consistent  random  inputs  can  then  be  treated  as  a  random  sample. 
The  hypotheses  tested  in  this  final  section  pertain  to  the  type  of 
probability  distribution  that  the  random  samples  most  closely  follow 
and  to  equality  of  population  variances  and  means. 


A.  Definition  of  the  Station-Keeping  Problem 

Even  for  an  accurate  nominal  trajectory,  unmodeled  as  well  as 
poorly  modeled  forces  on  the  spacecraft  will  generally  be  present,  and 
the  resulting  modeling  errors  may  be  a  contributing  cause  of  spacecraft 
drift  from  the  nominal  path.  As  was  discussed  in  Chapter  1,  a 
trajectory  near  a  collinear  libration  point  is  designed  to  excite  only 
the  stable  modes  associated  with  the  motion.  When  the  spacecraft 
deviates  from  this  nominal  trajectory,  the  unstable  mode  may  be 
excited,  and  drift  from  the  path  may  then  increase  exponentially  with 
time.  A  station-keeping  method  can  thus  be  used  to  combat  this  drift 
and  keep  the  spacecraft  acceptably  close  to  the  nominal  path.  Of 
course,  specifications  associated  with  the  station-keeping  procedure 
will  influence  the  deviations  from  the  nominal  path.  For  example,  as 
the  torus  about  the  nominal  path  is  relaxed  or  as  the  minimum 
acceptable  separation  time  between  control  inputs  is  increased,  drift 
from  the  nominal  path  may  increase. 

How  large  will  the  drift  be?  Table  2  contains  the  partial  results 
from  several  numerical  simulations  and  illustrates  the  magnitude  and 
direction  of  spacecraft  drift  from  a  nominal  halo  path  due  to 
deterministic  errors  in  initial  conditions  ("Ax")  or  solar  reflectivity 
("Ak" ).  For  these  simulations,  the  deterministic  initial  position 
errors  are  computed  from  a  consistent  position  on  the  halo  orbit.  The 
units  for  position  deviations  are  kilometers;  the  units  for  velocity 
deviations  are  in  millimeters  per  second. 
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Table  2.  Nominal  Path  Deviations  due  to  Deterministic  Errors. 

(The  units  for  position  are  kilometers;  velocity  is 
in  mm/sec. ) 


Deviation 

from  the 

Nominal  Path 

after 

Ax 

20  days 

40  days 

60  days 

80  days 

Ak=0 

Ak=. 13 

Ak=0 

Ak=. 13 

Ak=0 

Ak=. 13 

Ak=0 

Ak=. 13 

X 

1 

4.5 

14.6 

11.9 

59.0 

-5.9 

122.5 

-59.4 

242.5 

y 

1 

1.5 

-.6 

-.9 

-14.7 

-12.5 

-56.7 

-35.6 

-153.0 

z 

1 

2.1 

1.9 

2.1 

1.0 

5.6 

3.5 

6.9 

9.9 

X 

1 

3.3 

15.8 

6.8 

38.5 

-2.4 

64.3 

-28.4 

116.7 

y 

1 

-.5 

-3.8 

-1.9 

-13.0 

-.9 

-26.9 

5.0 

-61.5 

z 

1 

.  1 

-.2 

-.8 

-1.6 

-.2 

.1 

-.5 

6.4 

X 

1 

38.6 

138.2 

316.2 

677.0 

y 

1 

7.4 

-16.8 

-94.0 

-286.2 

z 

1 

15.0 

17.0 

13.8 

18.1 

X 

10 

44. 1 

83.9 

157.4 

317.1 

y 

10 

-4.2 

-24.6 

-58.9 

-151.5 

z 

10 

4.9 

-3.2 

-3.0 

10.4 

X 

1 

278.2 

930.2 

2253.8 

5025.0 

y 

1 

39.5 

-37.7 

-466.5 

-1616.8 

z 

1 

42.6 

177.8 

116.7 

100.5 

X 

100 

81.7 

538.0 

1088.6 

2324. 1 

y 

100 

-23.4 

-141.2 

-378.4 

-1051.2 

z 

100 

-18.5 

-18.9 

-43.5 

503.2 

X 

10 

20.2 

30.3 

51.7 

98.7 

85.6 

214.0 

141.7 

443.6 

y 

10 

7.3 

5.3 

-3.9 

-17.8 

-33.6 

-77.8 

-100.6 

-218.0 

z 

10 

7. 1 

6.9 

-1.9 

-3.  1 

-6.3 

-8.4 

-4.4 

-1.4 

X 

1 

11.6 

24. 1 

27. 1 

58.9 

39.8 

106.5 

63. 1 

208.3 

y 

1 

-4.0 

-7.3 

-9.0 

-20.0 

-16.1 

-42. 1 

35.8 

-102.3 

z 

1 

-4.3 

-4.5 

-6.3 

-7.2 

-3.1 

-2.8 

4.0 

11.0 
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Table  2,  continued. 


Ax 

I 

20  days 

leviatlon  from 
40  days 

the  Nominal  Path 
60  days 

after 

80  days 

Ak=0 

Ak=. 13 

Ak=0 

Ak=. 13 

Ak=0 

Ak=. 13 

Ak=0 

Ak=. 13 

X 

50 

90.0 

100.  1 

228.  1 

275.  1 

491.9 

620.3 

1035.4 

1337.4 

y 

50 

33.4 

31.5 

-17.5 

-31.3 

-127. 1 

-l,i.3 

-389.6 

-507.0 

z 

50 

29.6 

29.4 

-19.8 

-20.9 

-59.2 

-61.3 

-54.6 

-51.6 

X 

1 

48.4 

60.9 

117.5 

149.3 

227.4 

294. 1 

470.0 

615.3 

y 

1 

-19.7 

-23.0 

-40.3 

-51.4 

-83.5 

-109.5 

-217.2 

-283.7 

z 

1 

-23.4 

-23.7 

-31.0 

-31.8 

-16.0 

-15.7 

24.5 

31.4 

X 

50 

124.2 

354.3 

814.  1 

1772.0 

y 

50 

39.5 

-33.4 

-208 . 5 

-640. 1 

z 

50 

42.6 

-4.9 

-51.0 

-43.3 

X 

10 

81.7 

194.7 

387.2 

815.8 

y 

10 

-23.4 

-63.0 

-141.5 

-373.7 

z 

10 

-18.5 

-33.4 

-19.6 

35.4 

X 

100 

177.3 

187.4 

448.7 

495.8 

999.9 

1128.4 

2152.9 

2455.0 

y 

100 

66.2 

64.2 

-34.4 

-48.2 

-243.9 

-288. 1 

750.7 

-868.0 

z 

100 

57.7 

57.5 

-42.  1 

-43.3 

-125.3 

-127.4 

117.3 

-114.3 

X 

1 

94.5 

106.9 

230.5 

262.3 

462.0 

528.7 

1000.0 

1124.4 

y 

1 

-39.3 

-42.6 

-79.5 

-90.5 

-167.8 

-193.8 

-400.0 

-510.4 

z 

1 

-47.3 

-47.6 

-61.8 

-62.7 

-32.0 

-31.7 

-100.0 

57.  1 

X 

200 

351.7 

361.9 

890.0 

937.0 

2016.2 

2144.7 

4390.0 

4691.8 

y 

200 

131.6 

129.5 

-68.2 

-82.0 

-477.6 

-521.7 

-1472.3 

-1589.6 

z 

200 

1.1 

113.6 

-86.9 

-88.0 

-257.5 

-259.6 

-242.4 

-239.4 

X 

1 

186.0 

199.0 

456.6 

488.4 

931.4 

998.  1 

1998.6 

2144.0 

y 

1 

-78.6 

-81.9 

-157.8 

-168.9 

-336. 1 

-362. 1 

-897.0 

-963.4 

z 

1 

-95.  1 

-95.4 

-234.8 

-124.3 

-64.2 

-63.8 

-101.6 

1086.6 
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Notice  that  the  magnitude  of  the  drift  Is  sensitive  to  both 
initial  velocity  and  initial  position  deviations  and  that  a  solar 
reflectivity  constant  error  of  13%  also  has  a  substantial  effect  on  the 
drift.  Of  course,  the  drift  will  undoubtably  vary,  perhaps  in  a  minor 
way,  depending  on  the  size  and  shape  of  the  nominal  orbit  and  the 
location  on  the  path  at  which  the  deviations  are  evaluated.  Clearly, 
however,  a  spacecraft  can  quickly  move  quite  far  from  the  nominal  path. 

One  approach  that  might  be  considered  for  a  station-keeping 
strategy  is  to  redefine  the  nominal  path  such  that  it  passes  through 
the  current  position.  (Some  guidance  procedures  for  interplanetary 
missions  have  successfully  employed  such  a  scheme. )  Unfortunately,  an 
additional  complication  associated  with  a  libration  point  orbit  is 
that  a  nominal  path  cannot  generally  be  defined  through  all  possible 
positions.  A  bounded  orbit  may  not  exist  through  an  arbitrary  point, 
and  the  computational  difficulty  required  to  construct  a  nominal  orbit 
through  an  approximate  set  of  coordinates  also  makes  redefining  the 
nominal  path  during  flight  virtually  impossible.  Therefore,  a 
station-keeping  algorithm  that  returns  the  spacecraft  to  a  torus  about 
the  reference  trajectory  is  essential. 

A  station-keeping  strategy  must  combat  the  exhibited  drift  from 
the  nominal  path  while  satisfying  some  predetermined  set  of 
specifications.  It  is  not  uncommon  for  a  station-keeping  scheme  to 
include  constraints  related  to  quantities  such  as  the  timing  of 
manuevers,  control  magnitude,  and  deviation  distance  from  the  nominal 
path.  (Of  course,  additional  types  of  constraints  have  been 
implemented  in  other  guidance  schemes  and  could  possibly  prove 
beneficial  in  future  libration  point  orbit  studies. )  Identifying  a 
lower  bound  for  such  quantities  may  be  necessary  to  allow  time  for 
orbit  determination  and  to  help  ensure  efficient  use  of  control  energy. 
Figure  2-1  illustrates  the  decision  process  necessary  to  implement  the 
minimal  restrictions  typically  used  in  the  control  algorithms  developed 
in  this  chapter.  Each  of  the  decision  steps  depicted  in  Figure  2-1  may 
serve  to  delay  the  input  of  control  force. 
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For  a  collinear  libratlon  point  orbit,  a  small  deviation  from  the 
(unstable)  nominal  trajectory  can  lead  to  rather  large  drift  from  the 
path  in  a  short  time.  In  effect,  the  station-keeping  algorithm  must 
combat  both  the  current  drift  from  the  path  in  addition  to  the 
exponential  increase  in  the  drift  that  is  expected  if  no  correction  is 
implemented.  Any  delay  in  the  control  actuation  may  allow  the  drift  to 
increase  and  to  thus  compound  the  station-keeping  problem.  The 
magnitude  of  the  drift  can  be  clearly  seen  in  Table  2.  Of  course,  in 
such  a  nonlinear  problem,  it  is  possible  for  the  spacecraft  to  begin 
returning  to  the  nominal  path  on  its  own.  Another  possible  constraint 
for  the  station-keeping  scheme  could  be  a  check  on  the  growth  or  decay 
of  the  drift.  This  check  on  the  spacecraft’s  drift  from  the  nominal 
path  may  serve  to  delay  a  control  input,  or  it  may  prove  efficient  to 
input  control  energy  at  the  current  time  to  "assist"  the  drift  back  to 
the  path.  The  optimal  timing  (now  or  at  some  future  "best"  time)  to 
implement  a  given  control  input  is  a  possible  area  for  future  research. 

The  goal  of  the  station-keeping  routine  is  then  to  keep  the 
spacecraft  "close  enough"  to  the  reference  trajectory.  The  allowable 
deviations  may  depend  on  the  simulation  experience  with  a  given  control 
algorithm  and  on  mission  constraints,  including  the  propellant  cost 
that  can  be  tolerated.  The  evaluation  of  various  minimal  restrictions 
is  addressed  in  later  sections.  When  the  spacecraft  is  "near"  the 
nominal  trajectory,  it  is  reasonable  to  model  the  deviations  from  the 
reference  path  using  linear  analysis.  Consistent  with  such  a  model,  an 
investigation  of  the  problem  incorporating  linear  control  theory  is 
thus  Initiated. 


B.  Previously  Developed  Station-Keeping  Strategies 

Early  work  by  Dwivedi 125,261  resulted  in  development  of  a  control 
procedure  that  is  derived  from  minimization  of  a  cost  functional.  His 
approach  was  developed  for  interplanetary  space  flight;  however,  it 
also  shows  promise  as  a  controller  for  libratlon  point  orbits.  The 
cost  function  is  defined  by  weighting  both  the  control  energy  used  at 
an  initial  time  and  the  expected  deviations  of  the  spacecraft  from  the 
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nominal  path  at  two  distinct  times  in  the  future.  The  points  used  to 
compute  deviations  from  the  nominal  path  are  referred  to  in  this  work 
as  "target  points."  These  points  are  defined  along  the  nominal  path  at 
predetermined  times  that  are  called  "target  times."  (More  than  two 
target  points  can  be  used;  this  modification  proves  valuable  later  in 
this  eifort.j  A  version  uf  Dwivedi's  controller,  using  distinct 
manuever  times,  is  summarized  here,  and  an  extended  derivation  used  by 

[9j 

Pernicka  is  then  also  described. 

In  Dwivedi’s  algorithm,  the  cost  function  includes  several 
submatrices  partitioned  from  the  state  transition  matrix.  The  state 
transition  matrix  is  derived,  in  the  usual  way,  through  a  linearizing 
process  relative  to  the  nominal  path.  For  notatlonal  ease,  the  state 
transition  matrix  is  partioned  into  four  3x3  submatrices  as 


*(t  ,t  )  = 

k  0 


(2-1) 


Dwivedi’s  controller,  in  this  formulation,  computes  a  Av  input  (a 
3x1  vector),  with  magnitude  denoted  as  Av,  for  a  time  denoted  as  tQ. 
The  Av  is  added  to  the  initial  velocity  states  in  the  numerical 
integration  routine  in  order  to  change  the  deviation  of  the  spacecraft 
from  the  nominal  path  at  some  future  times.  In  this  derivation,  m^  is 
the  position  deviation  (a  3x1  vector)  and  v^  is  the  velocity  deviation 
(a  3x1  vector)  of  the  spacecraft  from  the  nominal  path  at  time  t  ,  with 
k  =  1,  2,  etc.  If  eQ  is  the  residual  velocity  (a  3x1  vector)  and  pQ  is 
the  residual  position  (a  3x1  vector)  relative  to  the  nominal  path  at 
time  t  ,  then  a  Av  input  at  t  could  be  used  to  predict  m^  for  k  =  1, 
2,  etc.  For  instance,  when  the  initial  position  xq  includes  an  initial 

velocity  perturbation  e  ,  a  delta  velocity  Av,  and  an  initial  position 

_  (24) 

perturbation  pQ,  the  state  propagation  equation  provides: 
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(2-2) 


x 

k 


m 

p 

k 

=  *(t  t  )  X  =  *(t  ,t  ) 

0 

V 

K  U  U  K  U 

e  +Av 

k 

0 

As  an  example,  for  two  target  times  t  and  t2,  and  by  assuming  that  pQ 
is  the  zero  vector, 


m  =  B  e 
1  10  o 


B  Av 
10 


(2-3) 


and 


m  =  B  e  +  B  Av. 
2  20  o  20 


(2-4) 


Dwivedi’s  approach  was  developed  for  application  to  interplanetary 
missions  and,  consequently,  includes  the  assumption  that  position 
deviations  from  the  nominal  path  have  a  minor  effect  on  spacecraft 
drift  when  compared  to  velocity  deviations.  Hence,  pQ  is  allowed  to  be 
the  zero  vector  when  deriving  equations  (2-3)  and  (2-4)  from  equation 
(2-2).  The  target  times  t^  and  t  are  computed  by  adding  incremental 
times  At^  and  At2  to  tQ,  respectively.  For  instance,  if  At^  =  40  days, 
At  =  70  days,  and  t  =  0  days,  then  t  =  t  +  At  =40  days  and  the 
second  target  time  t  =  t  +  At2  =  70  days. 

The  cost  function  that  will  be  minimized  is 


J(Av)  =  AvT  Q  Av  +  mT  R  m  +  mT  S  m  (2-5) 

112  2 

where  the  weighting  matrix  Q  is  symmetric  positive  definite,  and  the 
weighting  matrices  R  and  S  are  positive  semidef inite.  Equation  (2-5) 
can  be  written  in  terms  Av  by  substituting  equations  (2-3)  and  (2-4) 
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into  equation  (2-5).  Figure  2-2  depicts  a  version  of  Dwivedi’s 

station-keeping  routine  that  allows  the  control  input  at  a  specified 

time  other  than  t  and  that  uses  two  target  times,  which  are  actually 

incremental  times  At  and  At  added  to  the  initial  time  t  . 

12  o 

Determination  of  the  Av  corresponding  to  the  relative  minimum  of 
this  coat  function  allows  a  linear  equation  for  the  optimal  control 
input  (Av  )  to  be  found:  125,261 


Av 


-[Q  +  BT  R  B  +  BT  S  B  ]_1[BT  R  B  +  bt  S  B  ) 
10  10  20  20  10  10  20  20 


(2-6) 


Note  that  equation  (2-6)  assumes  control  implementation  at  time  t  . 
This  derivation  could  be  generalized  to  include  the  possibility  of 
time-varying  weighting  matrices  and  a  Av  input  at  some  time  t  after  t  , 
and  these  generalizations  could  be  the  subjects  of  valuable  future 
research  concerning  libratlon  point  orbital  control. 

1 9  27] 

Howell  and  Pernicka  ’  modified  the  above  controller  to  also 
include  the  effects  of  position  deviations  at  time  t  .  Dwivedi’s 
approach  assumes  that  velocity  perturbations  have  a  much  greater 
propagative  effect  than  position  deviations  from  the  nominal  path. 
Because  of  the  unstable  nature  of  libration  point  orbits,  it  was 
reasoned  that  both  position  and  velocity  residuals  from  the  nominal 
path  should  be  included  in  the  error  propagation  equations  such  as 
equation  (2-3)  and  equation  (2-4)  for  Dwivedi’s  controller.  An  example 
of  the  relationship  between  position  and  velocity  errors  at  t  and  the 
resulting  residuals  at  varying  time  steps  later  was  illustrated  in 
Table  2. 


40 


*1  *2 

m^S  m1  rn^S  in, 

e  =  VELOCITY  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t.  THIS  IS  A  3  x  1  VECTOR. 

POSITION  DEVIATIONFROM  THE  NOMINAL 

PATH  AT  TIME  IF  A  v  IS  INPUT  AT  TIME  t. 
THIS  IS  A  3  x  1  VECTOR. 

m^=  POSITION  DEVIATIONFROM  THE  NOMINAL 
PATH  AT  TIME  t2  IFav  IS  INPUT  AT  TIME  t. 
THIS  IS  A  3  x  1  VECTOR. 

R,  S  ARE  3x3  POSITIVE  SEMIDEFINITE 
WEIGHTING  MATRICES  USED  IN  THE  COST 
FUNCTION. 


Figure  2-2.  Dwivedi  Control  Routine. 
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Certainly,  both  position  and  velocity  offsets  from  the  nominal 
trajectory  do  appear  to  make  a  substantial  contribution  to  spacecraft 
drift  from  the  nominal  path.  For  this  derivation,  pQ  is  the  nonzero 
position  residual  vector  (3x1)  and  eQ  is  again  the  nonzero  velocity 
residual  vector  (3x1)  at  time  tQ;  by  using  equation  (2-2),  the 
following  two  equations  can  be  derived: 


m  = 

l 


B  e  + 
10  o 


B  Av  +  A  p 
10  10*0 


(2-7) 


B  e 
20  0 


B  Av 
20 


+  A  p  . 
20  *0 


(2-8) 


Figure  2-3  is  a  depiction  of  the  Howell/Pernicka  controller  that  uses  a 
control  input  at  t  and  also  uses  two  target  times  (at  predetermined 
incremental  times  beyond  t  ).  The  cost  function  that  is  minimized 
remains  as 

J(Av)  =  AvT  Q  Av  +  mT  R  m  +  mT  S  in  .  <.2-9) 

112  2 

Using  equations  (2-7)  and  (2-8)  in  equation  (2-9),  the  optimal  control 
is  obtained  by  minimizing  the  cost  function  in  terms  of  Av,  and  it  then 

.  [9,27] 

becomes 


Av*=  -  [  Q  +  BT  R  B  +  BT  S  B  ]'l[(BT  R  B  +  BT  S  B  )  e  + 
10  10  20  20  10  10  20  20  0 


(BT  R  A  +  BT  S  A  )  p  )  .  (2-10) 

10  10  20  20  0 
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pi  cpl 


rr^Rm,  r?£S  rr^ 


e0=  VELOCITY  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t0 .  THIS  IS  A  3  x  1  VECTOR. 

p0=  POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t0 .  THIS  IS  A  3  x  1  VECTOR. 

POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  IF  Av  IS  INPUT  AT  TIME  V 
THIS  IS  A  3  x  1  VECTOR. 

m  =  POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t2  IFAv  IS  INPUT  AT  TIME  ^ 
THIS  IS  A3  xl  VECTOR. 


R,  S  ARE  3x3  POSITIVE  SEMIDEFINITE 
WEIGHTING  MATRICES  USED  IN  THE  COST 
FUNCTION. 


Figure  2-3.  Howell/Pernicka  Control  Scheme 
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[9  27] 

Clearly,  the  Howell/Pernicka  formula  ’  for  the  optimal  control 

input  is  slightly  more  complicated  than  that  derived  using  Dwivedi’s 

approach.  The  inclusion  of  the  position  deviation  (pQ)  in  equation 

(2-10),  versus  using  equation  (2-6),  may  be  shown,  in  fact,  to  reduce 

the  control  costs.  Both  the  Dwivedi  and  the  Howell/Pernicka  control 

algorithms  can  be  used  for  quasi-perlodic  and  periodic  libration  point 

[9  27] 

orbits.  Howell  and  Pernicka  ’  used  their  control  law  (2-10)  for  a 
spacecraft  in  a  halo-type  orbit  near  the  interior  libration  point  in 
the  Sun-Ear th+Moon  system,  and  their  algorithm  was  formulated  to 
include  several  features.  Minimal  separation  times  of  up  to  80  days 
between  control  inputs  were  used.  This  is  a  realistic  feature  because 
the  orbit  determination  process  and  the  control  input  computations 
require  some  minimal  work  time.  The  80-day  control  separation  time  was 
selected  to  roughly  correspond  to  the  timing  between  manuevers  used  for 
ISEE-3.  They  also  included  a  minimal  control  input  ("Av")  magnitude  in 
their  formulation.  Modern  propellant  devices  do  have  restrictions 
concerning  the  minimal  control  energy  that  can  be  accurately  expended, 
and  the  errors  relative  to  the  commanded  control  may  be  a  function  of 
the  control  magnitude.  (In  fact,  the  control  uncertainty  modeled  later 
in  the  station-keeping  simulations  may  be  inversely  related  to  control 
magnitude. )  The  values  that  Howell  and  Pernicka  used  for  the  minimal 
Av  ranged  between  .01  and  .5  meters  per  second  depending  on  the 
specific  simulation  completed. 

Incorporating  a  torus  of  acceptable  size  about  the  nominal  path  is 
also  a  useful  feature.  Howell  and  Pernicka  used  a  torus  that  ranged 
between  0  and  100  kilometers  in  radius.  If  the  specific  torus 

dimensions  were  not  violated,  even  when  all  other  conditions  for  a 

control  input  were  met,  no  Av  would  be  expended.  The  control  costs 

(9  27) 

found  in  Howell  and  Pernicka  ’  are  generally  comparable  to  those 

found  in  other  investigations.  (These  results  will  be  included  in  the 
last  section  of  this  chapter. )  However,  it  is  noted  that  simulations 
that  included  tracking  and  injection  errors  and  a  minimum  control 
separation  time  of  80  days  had  relatively  high  propellant  costs.  These 
propellant  costs  are  also  obviously  a  function  of  the  target  point 

spacing;  therefore,  changes  in  specifications  may  enable  some  reduction 
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in  the  cost.  It  should  also  be  noted  that  ISEE-3  was  successfully 
controlled  at  a  cost  that  was  lower  than  the  propellant  cost  predicted 
by  pre-mission  station-keeping  simulations.  Further  investigation  of 
the  choices  of  weighting  matrices  {both  constant  ana  time-varying),  the 
number  of  target  points,  and  the  target  spacing  may  prove  valuable; 
alternatively,  a  controller  that  also  weights  velocity  residuals  at  the 
target  times  may  prove  to  be  an  Improvement. 


C.  Delta-Velocity  Controller  I 

An  innovation  that  also  adds  velocity  residual  weighting  in  the 

[9  27] 

cost  function  of  the  Howel 1/Pernicka  ’  controller  provides 

promising  results.  For  this  controller,  the  velocity  residuals  at  t 
(denoted  by  v^ )  and  t^  (denoted  by  v^)  are  3x1  vectors  that  can  be 
computed  as 


(2-11) 


The  cost  function  for  this  controller  (Delta-Velocity  Controller  I)  is 
then 

J(Av)  =  AvT  Q  Av  +  mT  R  m  +  vT  R  v  +  mT  S  m  +  vT  S  v  ,  (2-13) 

1  1  1  v  1  2  2  2  V  2 
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where  the  added  weighting  matrices,  R  and  S  ,  are  positive 

V  V 

semidef inite.  (This  formulation  is  further  complicated  by  the  addition 

of  two  weighting  matrices  that  must  be  somehow  chosen.  In  general, 

there  does  not  appear  to  be  an  strategy  that  can  optimally  select  the 

weighting  matrix  entries  for  this  controller  or  for  the  station-keeping 

routines  of  Dwivedi  or  Howell  and  Pernicka.  Extensive  experimentation 

in  addition  to  investigator  judgement  are  often  successfully  used; 

however,  a  selection  method  for  the  weighting  matrix  entries  will  be  a 

vali. b,e  area  of  future  research.) 

Using  equations  (2-11)  and  (2-12),  the  following  equations  for  the 

residuals  at  time  t  and  t  can  be  found: 

1  2 


m=B  e  +  B  Av  + 
1  10  o  10 


A  p  , 
10  o 


(2-14) 


m  =  B  e  +  B  Av  + 
2  20  o  20 


A  p  , 
20  o 


(2-15) 


v  =  D  e  +  D  Av  + 
1  10  o  10 


C  p  , 
10r0 


(2-16) 


v  =  D  e  +  D  Av  +  C  p.  (2-17) 

2  20  0  20  20  o 


Recall  that  t^  and  t  are  computed  by  adding  incremental  times  Atj 

and  At^,  respectively,  to  the  initial  time  tQ.  Figure  2-4  is  a 

depiction  of  the  Delta-Velocity  Controller  I  that  uses  a  control  input 

at  the  initial  time  t  and  Includes  two  target  times  (t  and  t  ).  The 

o  12 

weighting  of  velocity  errors  in  the  cost  function  will  obviously 
complicate  the  derivation. 
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loP  |d? 


^  2 

rr£S  rr^ 
Vg  Svv2 

e"0=  VELOCITY  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  ^  THIS  IS  A  3  x  1  VECTOR. 

p0=  POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t0 .  THIS  IS  A  3  x  1  VECTOR. 

m,  =  POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t!  IF  AV  IS  INPUT  AT  TIME  t0. 

THIS  IS  A  3  x  1  VECTOR. 

vj  =  VELOCITY  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t,  IF^  IS  INPUT  AT  TIME  t0- 
THIS  IS  A  3  x  1  VECTOR. 

mp=  POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t2  IF AV  IS  INPUT  AT  TIME  t0  • 
THIS  IS  A  3  x  1  VECTOR. 

v2=  VELOCITY  DEVIATION  FROM  THE  NOMINAL 

PATH  AT  TIME  t2IFAV  IS  INPUT  AT  TIME  t0- 
THIS  IS  A  3  x  1  VECTOR. 

R,  S,  Rv,  SVARE  3x3  POSITIVE  SEMIDEFINITE 
WEIGHTING  MATRICES  USED  IN  THE  COST 
FUNCTION 

Figure  2-4.  Delta- Velocity  Controller  I. 


m[R  m, 

v;  rvvi 
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Equations  (2-14)  through  (2-17)  can  be  substituted  into  equation 
(2-13),  and  the  relative  minimum  of  this  function  of  Av  (after 
considerable  algebra)  is  found  to  be 


Av  =  -  [  Q  + 


BT  R  B  +  BT  S  B  +DT  R  D  +  DT  S  D  ]_1x 

10  10  20  20  10  v  10  20  v  20 


[  (BT  R  B  +  BT  S  B  +  DT  R  D  +D  S  DT  )  e  + 

10  10  20  20  10  v  10  20  v  20  0 


(BT  R  A  +  8T  S  a  +  DT  R  C  +DT  S  C  )  p  ] .  (2-18) 

10  10  20  20  10  v  10  20  v  20  *0 


In  the  final  section  of  the  chapter,  this  control  law  will  be 
shown  to  provide  excellent  results  for  station-keeping  of  a  libration 
point  orbit.  However,  when  the  minimum  separation  time  between  control 
inputs  is  increased  to  60  or  80  days,  this  controller  tends  to  exhibit 
an  increase  in  cost.  Perhaps  a  controller  that  looks  further  downtrack 
may  provide  some  improvement  in  control  costs. 


D.  Delta-Velocity  Controller  II 

One  way  to  add  cost  function  weighting  to  residuals  farther  along 

the  track  would  be  to  increase  the  size  of  incremental  times  At  and 

l 

At^  Alternatively,  a  third  target  point  can  be  added  to  the 
formulation  of  Delta-Velocity  Controller  I;  this  adjustment  can  permit 
the  target  times  to  be  approximately  equally  spaced  (an  arbitrary 

choice)  with  A^  =  40  days,  for  instance.  This  adjustment  will  be 
shown  to  increase  the  robustness  and  decrease  the  cost  of  the 

station-keeping  algorithm,  especially  when  the  minimum  control  input 
separation  times  are  extended  to  60  or  80  days. 
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For  this  formulation,  position  and  velocity  residuals  at  time  t_ 

are  added  as  3x1  vectors  m  and  v  to  give: 

3  3 


m 

l 

A 

10 

B 

10 

P0 

v 

l 

C 

10 

D 

10 

e  +  Av 

0 

(2-19) 


m 


A  B 
20  20 


C  D 
20  20 


e  +  Av 
o 


(2-20) 


in 

3 

A 

30 

B 

30 

P0 

V 

3 

c 

30 

D 

30 

e  +  Av 

0 

(2-21) 


Figure  2-5  depicts  the  elements  of  Delta-Velocity  Controller  II  that 
includes  the  weighting  of  velocity  errors  at  the  target  points,  the  use 
of  three  target  times,  and  a  control  input  at  time  t  . 

The  cost  function  now  Includes  terms  for  the  third  time  step  and 
becomes 


J(Av)=AvTQ  Av  +  mTR  in  +  vTR  v  +  mTS  in  +  vTS  v  +  mTT  m  +  vTT  v 

1  1  lvl  2  2  2V2  3  3  3v 


v‘T  v  ,  (2-22) 

3  v  3 


49 


laf  \& 


mfRr^  nf^  mjT 

v/  RvVj  v2t  Sv  v2  v3t  Tv 


e^=  VELOCITY  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  ^  .  THIS  IS  A  3  x  1  VECTOR. 

pU  POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t0.  THIS  IS  A  3  x  1  VECTOR. 

m=  POSITION  DEVIATION  FROM  THE  NOMINAL 
1  PATH  AT  TIME  ttIFAv  IS  INPUT  AT  TIME  to- 
THIS  IS  A  3  x  1  VECTOR. 

7=  VELOCITY  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t,  IFav  IS  INPUT  AT  TIME  to- 
THIS  IS  A  3  x  1  VECTOR. 

m2=  POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t2  IF  3?  IS  INPUT  AT  TIME  tQ 
THIS  IS  A  3  x  1  VECTOR. 

~2  =  VELOCITY  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t,  IFAv  i  INPUT  AT  TIME  to 
THIS  IS  A  3  x  1  VECTOR. 

m3=  POSITION  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  t3  IF  3?  IS  INPUT  AT  TIME  tQ. 
THIS  IS  A  3  x  1  VECTOR. 

V3=  VELOCITY  DEVIATION  FROM  THE  NOMINAL 
PATH  AT  TIME  ^  IF  ^  IS  INPUT  AT  TIME  V 
THIS  IS  A  3  x  1  VECTOR. 

R,  Rv ,  S,  Sv  ,  T,  \  ARE  3x3  POSITIVE  SEMI- 

DEFINITE  COST  FUNCTION  WEIGHTING  MATRICES 


Figure  2-5.  Delta-Velocity  Controller  II. 
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J3' 


where  the  added  weighting  matrices,  T  and  T  ,  are  positive 

semidef inite.  The  cost  function  can  be  written  in  terms  of  Av  by  using 

substitutions  for  m  and  v  ,  with  k  =  1,  2,  and  3,  derived  from 

k  k 

equations  (2-19),  (2-20),  and  (2-21).  The  relative  minimum  of  this 

function  of  Av  (after  considerable  algebra)  is  found  to  be 


Av*=-[Q+BT  r  b  +bt  s  b  +bt  t  b  +dt  r  d  +dt  s  d  +dt  t  d 

10  10  20  20  30  30  10  v  10  20  v  20  30  v  30 


-1 

X 


( 


(BT  R  B  +bt  S  B  +bt  T  B  +dt  R  D  +dt  S  D  +dt  T  D 

10  10  20  20  30  30  10  v  10  20  v  20  30  v  30 


+ 


(bt  r  a  +bt  s  a  +bt  t  a  +dt  r  c  -*-dt  s  c  +dt  T  C  )p  ] . 

10  10  20  20  30  30  10  v  10  20  v  20  30  v  30  0 


(2-23) 


The  general  method  described  above  could  certainly  accommodate  the 
inclusion  of  additional  target  points.  The  relative  value  of  using 
additional  target  points  and  an  algorithm  for  selecting  the  weighting 
matrix  entries  both  seem  to  be  potential  areas  for  future  research. 
Preliminary  work  using  entries  from  the  state  transition  matrices  has 
shown  some  degree  of  promise  in  choosing  the  weighting  matrix  entries. 
Results  to  date  are  inconclusive. 

The  third  type  of  controller  investigated  is  termed  an  on/off 
station-keeping  method  because  the  control  energy  is  input  as  an 
acceleration,  added  for  a  few  days  and  then  removed  for  varying  periods 
up  to  80  days.  The  on/off  type  of  controller  is  fundamentally  a 
modification  of  a  discrete-time  continuous  controller  that  is  modified 
to  be  more  operationally  feasible. 
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E.  On/Off  Controller 


The  station-keeping  method  considered  for  use  here  is  state 
feedback  control  computed  from  minimization  of  a  quadratic  cost 

[28-31 J  [32  33] 

function.  Other  researchers  ’  have  used  a  similar  control 
scheme  for  work  concerning  planned  libration  point  orbits  in  the 
Earth-Moon  and  Sun-Earth  systems.  The  method  has  been  shown  to  produce 
competitively  low  control  costs;  however,  this  type  of  scheme  has 
significant  drawbacks  for  actual  implementation.  The  method  as  used 
here  assumes  piecewise  constant  control  inputs,  yet  thrusters  typically 
are  not  designed  for  constant  operation.  (This  type  of  thruster  is 
being  developed  now,  however,  for  use  in  the  next  decade. )  Generally, 
impulsive  control  inputs  are  preferred  in  practice. 

The  control  strategy  considered  here  incorporates  the  use  of  a 
torus  about  the  reference  trajectory;  the  control  input  is  delayed 
until  the  limits  of  the  predetermined  torus  are  violated.  Because  of 
actual  mission  constraints,  this  formulation  also  includes  a  specified 
minimum  number  of  days  between  manuevers  and  a  minimum  control  input 
magnitude.  Simplifying  assumptions,  related  to  accommodation  of  the 
time-varying  nature  of  the  residuals,  will  be  discussed  later.  All  of 
these  minimal  constraints  and  as  yet  unspecified  simplifying 
assumptions  make  this  formulation  truly  "suboptimal"  yet 
computationally  simple;  however,  a  great  deal  of  problem  insight  can  be 
gained  from  this  analysis. 

Throughout  the  formulation  of  all  of  the  station-keeping 
algorithms,  the  reference  trajectory  is  defined  by  six  cubic  splines, 
one  for  each  of  the  six  states.  The  continuous,  linearized  residual 
model  is  given  in  the  following  form: 


x ( t )  =  A C t )  x(t)  +  B  u(t) 


(2-24) 


u(t)  =  G(t)  x(t), 


(2-25) 
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where 


B  = 


0  0  0 
0  0  0 
0  0  0 
1  0  0 
0  1  0 
0  0  1 


(2-26) 


and  A(t)  is  the  6x6  Jacobian  matrix  derived  in  Chapter  1.  The  state 
feedback  gain  matrix  is  G(t).  For  computation  of  the  piecewise 
constant  suboptimal  gains,  the  system  (2-24)  is  discretized  to 

j  [29-31,34] 

produce: 


x  =  *(k+l,k)  x  +  B  u(k), 

k+l  k  D 


(2-27) 


where  x^  is  the  residual  state  vector  at  time  step  k;  *(k+l,k)  is  the 
state  transition  matrix  at  time  step  k+l  relative  to  time  step  k;  B^  is 
the  discrete  matrix  derived  from  B;  and  u(k)  =  G  x  is  the 

k  k 

piecewise  constant  state  feedback  control.  The  control  energy  is 
computed  in  an  optimal  way  through  minimization  of  a  quadratic  cost 
function. 

The  feedback  gain,  G^,  can  be  calculated  by  minimizing  the 
following  total  cost  function: 


n- 1 

V(p)  -  E  [xj  p  Q  X  ♦  uT(k)  R  u(k)], 

k  =  0 


(2-28) 
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where  Q  is  the  diagonal  weighting  matrix  for  state  residuals;  R  is  the 
diagonal  weighting  matrix  for  control  input;  n  is  the  number  of  time 
steps  used  for  the  control  inputs;  p  is  the  scalar  that  can  be  used  to 
vary  "tightness";  and  terminal  state  constraints  are  omitted. 

The  number  of  time  steps  needed  for  a  given  trajectory  is 
determined  by  the  time  in  days  between  each  computed  change  in  control 
effort  and  by  the  duration  of  the  orbit.  For  example,  a  6-year 
trajectory  with  control  inputs  computed  in  20-day  segments  would 
require  n  =  108  control  time  steps.  In  general,  minimization  of  the 
cost  function  V(p)  is  obtained  through  a  sequence  of  difference 

[30  31] 

equations  that  can  be  solved  by  a  backward  sweep  method.  ’ 

One  further  computational  simplification  employed  here  is  to  use 
the  stabilizing  steady-state  gain  solution  to  the  cost  function  V(p). 
This  steady-state  assumption  (n-x»)  yields  the  following  matrix  solution 

[29] 

for  the  state  feedback  gain  matrix  G^: 


G  =  -(R  +  BT  K  B  rV  K  *(k+l,k), 

k  0  k  D  D  k 


(2-29) 


where 


K  =$T (k+1 ,  k)K  4>(k+l  ,k)- 

k  k 


$T  (k+1  ,k)K  B  (R+BT  K  B  TV 

k  D  D  k  D  D 


K  #(k+l ,k)+pQ, 

k 


(2-30) 


u 

k 


=  G  x  =  (u  , 

k  k  1 


(2-31) 
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Equation  (2-30)  is  the  algebraic  matrix  “Riccati"  equation  used  to 

compute  optimal  steady-state  feedback  gains.  This  algebraic  matrix 

"Riccati"  equation  is  a  simplification  of  a  matrix  difference  equation 

that  has  a  continuous  counterpart — the  "Riccati"  matrix  differential 

equation  used  to  compute  optimal  state  feedback  in  continuous 

time-varying  systems.  These  matrix  "Riccati"  equations  were  so  named 

1 35] 

by  Rudolf  E.  Kalman  in  1960.  In  1724,  Jacopo  Francesco,  Count 

Riccati  (1676-1754)  of  Venice,  had  considered  the  solution  of  a  special 

[  36  3*7  ] 

form  of  a  scalar  differential  equation  in  his  work  in  acoustics. 

His  equations  appeared  in  many  applications  related  to  Bessel 
functions.  Some  20  years  earlier,  James  Bernoulli  (1654-1705)  worked 

[37] 

on  solutions  to  a  similar  differential  equation.  Jean  Le  Rond 

D’Alembert  (1717-1783)  was  the  first  to  work  on  a  more  general  form  of 

these  scalar  differential  equations  and  used  the  name  "Riccati 

[36  ] 

equation"  for  the  equations  of  that  general  form. 

In  this  work,  the  continuous  controller  uses  piecewise  constant 
gains.  These  steady-state  gains  are  computed  for  time  steps  of  20  days 

[34  381 

in  a  suboptimal  scheme.  ’  ‘  The  exact  equations  of  motion  are 

integrated,  Incorporating  the  control  inputs,  in  the  elliptic 

restricted  three-body  model  using  a  Runge-Kutta  fifth-order  numerical 
integration  routine. 

The  differential  equations  to  be  integrated  are  then 


x  -  2  0  y  =  U  +  9y  +  u, 

X  1 


(2-32) 


y  +  2  0  x  =  U 

y 


(2-33) 


z  =  U  +  u  . 

z  3 


(2-34) 


55 


The  control  energy  is  input  through  the  constant  accelerations 
denoted  in  equations  (2-32),  (2-33),  and  (2-34)  as  u  ,  u  ,  and  u  , 

12  3 

respectively.  This  state  feedback  contr '  ler,  using  steady-state 
gains,  could  be  formulated  to  compute  a  new  control  u^  every  control 
time  step.  The  control  time  step  might  be  every  20,  30,  45,  or  even  60 
days.  In  this  way,  a  low  level  of  accelerative  control  energy, 
changing  periodically,  could  be  used  continuously  to  maintain  the 
spacecraft  within  a  torus  about  the  nominal  path.  Figure  2-6  depicts 
the  control  scheme  for  the  discrete-time  continuous  controller.  This 
type  of  continuous  controller  formulation  has  shown  excellent  results; 
however,  the  very  low  level  of  commanded  control  input  is  not 
operationally  feasible. 

A  more  practical  station-keeping  algorithm  could  use  an  on/off 
control  scheme  that  would  also  incorporate  the  various  minimal  control 
constraints  mentioned  previously.  Even  though  an  impulsive 

(delta-velocity)  controller  is  generally  considered  to  be  the 
currently  preferred  method  for  station-keeping,  investigation  of  an 
on/off  conti  oiler  can  lead  to  valuable  problem  insight  and  may  someday 
prove  useful,  given  future  technological  advances. 

For  the  On/Off  Controller,  the  control  inputs  are  set  at  a 

constant  magnitude  for  a  given  20-day  time  period  and  are  then  off 

(zero  magnitude)  until  certain  minimal  constrants  are  met  or  exceeded. 

In  this  formulation,  the  control  energy  is  zero  unless  a  specified 

acceptable  deviation  distance  from  the  nominal  path  is  violated,  a 

minimum  control  separation  time  is  exceeded,  and  a  minimum  control 

magnitude  is  surpassed.  The  On-Off  Controller  is  depicted  in  Figure 

2-7.  The  control  effort  is  computed  by  using  the  spacecraft’s  position 

and  velocity  errors  relative  to  the  nominal  state  at  time  t^  to 

calculate  u  =G  x  where  G  is  the  constant  state  feedback  gain  matrix 
k  k  k  k 

computed  in  equation  (2-29)  for  the  given  control  time  step,  and  x^  is 
the  state  vector  of  residuals  from  the  nominal  path. 
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to  t  j  t2  t3  t4  t  5 

7,  =[u,,u2,u3  ] . ,  i  =  1 ,  2,  3 . 

|TT|  =  [u?+u|  +  u|]’'2 

1 1  '  to  =  ^  2  ~  ^1  =  ^3"  ^2=^4  "  ^3  =*  5  "  ^  4  =  A  t 

At  COULD  BE  20  DAYS  OR  30  DAYS  OR  40  DAYS 
FOR  THIS  RESEARCH,  At  IS  20  DAYS. 


Figure  2-6.  Discrete  Time  Continuous  Controller. 


57 


At  COULD  BE  20  DAYS  OR  30  DAYS  OR  40  DAYS 


FOR  THIS  RESEARCH,  At  IS  20  DAYS. 


Figure  2-7.  On/Off  Controller. 


58 


The  weighting  factors  (p,  Q,  and  R)  in  equation  (2-28)  determine 
the  contribution  to  the  total  cost  of  deviations  from  the  nominal  path 
as  compared  to  the  costs  of  the  control  inputs  required  to  counter 
these  deviations.  The  positive  scalar  weight,  p,  has  been  used  in 

129  31  33] 

other  research  ’  ’  and,  for  this  effort,  p  determines  the 
“tightness"  of  the  control;  that  is,  a  relatively  large  value  of  p 
causes  residuals  relative  to  the  nominal  path  to  be  more  costly.  This 
will,  in  general,  force  the  control  algorithm  to  input  comparatively 
larger  control  effort  to  maintain  the  spacecraft  closer  to  the  nominal 
path.  The  R  weighting  matrix  in  the  cost  function,  V(p),  in  (2-28)  is 
arbitrarily  chosen  as  the  3x3  identity  matrix.  This  value  for  R  gives 
the  control  input  in  each  direction  equal  weighting. 

The  elements  in  the  6x6  diagonal  matrix  Q  may  be  chosen  in  several 
ways.  One  particular  method  for  choosing  these  entries  provided 
results  that  were  inconclusive  here  yet  may  prove  valuable  in  future 
research.  For  each  control  time  step  (20  days),  the  entries  in  the  Q 
weighting  matrix  were  computed  in  a  predetermined  way  from  elements  of 
the  state  transition  matrix  evaluated  at  that  time  step.  Use  of  the 
sensitivities  reflected  in  the  state  transition  matrix  slightly 
decreased  the  total  control  cost  corresponding  to  this  station-keeping 
method.  This  selection  method  for  the  elements  of  the  Q  weighting 
matrix  is  an  area  for  future  investigation  and  is  not  presented  here. 

Secondly,  both  the  continuous  and  the  on/off  control  methods  can 
incorporate  the  estimated  off-course  deviation  from  the  nominal  path  at 
the  end  of  a  control  time  step  in  order  to  modify  the  control  effort 
input  at  the  start  of  the  step.  This  "shooting  method"  uses  the 
otherwise  planned  control  input  to  predict  the  resulting  error  at  the 
end  of  the  control  time  step  or  at  some  point  further  along  the  track, 
and  then  adjusts  the  initial  gains  to  accommodate  the  predicted  errors. 
This  modification  has  been  shown  to  substantially  decrease  the  total 
control  costs  as  compared  to  other  adjustments  investigated  in  this 
research.  The  results  of  several  differing  station-keeping  methods  are 
now  presented  in  Chapter  3. 


59 


CHAPTER  3:  STATION-KEEPING  RESULTS  AND  COMPARISONS 


Ideally,  the  control  scheme  modifications  that  art  investigated 
here  and  the  overall  results  from  this  preliminary  work  should  be 
compared  to  other  Lagrange  point  station-keeping  results.  A  criteria 
for  this  comparison  could  be  the  magnitude  of  the  total  control  effort 
that  is  needed  to  successfully  maintain  the  vehicle,  satisfying  all 
station-keeping  requirements,  for  a  1 ’.brat ion  point  orbit  of  some  fixed 
duration.  The  fixed  duration  for  these  comparisons  will  be  a  2-year 
segment  of  the  nominal  trajectory.  The  results  of  the  various 
station-keeping  approaches  that  are  evaluated  in  this  research  will  be 
shown  to  produce  excellent  results.  However,  before  presenting  any 
results  and  then  discussing  the  comparisons,  it  is  necessary  to  address 
two  additional  questions: 

1.  Are  the  Av  random  variables  (approximately)  normally 
distributed? 

2.  What  is  a  reasonable  sample  size? 

In  this  section,  the  results  of  each  2-year  station-keeping 
simulation  is  considered  a  function  of  several  random  inputs  and  is, 
therefore,  treated  as  a  random  trial.  Each  station-keeping  simulation 
produces  a  total  propellant  value  (Av^)  that  is  a  scalar  measure  of  the 
cost  of  station-keeping  for  that  2-year  period.  The  simulation  that 
produces  the  A*'  variable  is  subject  to  several  random  inputs  and, 
consequently,  will  vary  for  each  Monte  Carlo  trial  (simulation)  of  the 
station-keeping  algoritm.  Several  independent  simulations  (random 
trials),  where  the  input  random  variables  have  consistent  statistical 
characteristics  (mean  and  variance),  will  produce  a  random  sample  of 
Av^  results  (one  for  each  simulation).  That  is,  each  station-keeping 
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run  is  treated  as  a  random  trial  with  random  variable  Av  (the 

T 

magnitude  of  total  Av  for  the  two-year  orbit),  and  a  group  of  random 
trials  for  which  all  the  input  variables  (torus  size,  tracking  error 
levels,  etc)  remain  the  same  are  then  treated  as  a  random  sample.  This 
resulting  sample  can  then  be  tested  to  see  if  it  fits  the  Gaussian 
distribution;  the  somewhat  related  question  of  the  required  number  of 
station-keeping  runs  per  sample  must  then  be  addressed. 

If  the  results  of  the  station-keeping  simulations  were  simply  a 
linear  combination  of  Gaussian  random  variables,  the  resulting 
distribution  of  Av^  would  be  Gaussian.  This  simple  representation  is 
not  true  here;  therefore,  it  is  necessary  to  complete  a  chi-squared 
goodness  of  fit  test  on  the  results  of  several  trial  runs  in  an  initial 
sample.  (Certainly,  there  are  other  statistical  tests  that  can  be 
used  to  determine  whether  a  sample  is  drawn  from  a  Gaussian 
distribution;  the  chi-squared  test  is  not  the  most  powerful,  but  it  is 
the  most  easily  presented. )  Choosing  a  sufficiently  large  sample  size, 
comparing  population  means,  and  constructing  confidence  intervals  can 
be  simplified  by  knowing  the  probability  distribution  of  the  propellant 
costs.  Hence,  the  computation  of  a  sufficient  sample  size  apriori  will 
permit  the  comparison  of  station-keeping  methods  with  a  minimal  number 
of  simulation  runs. 

The  results  and  comparisons  that  follow  also  include  statistical 
comparison  tests  for  equivalent  variances  and  equivalent  means  of  the 
values  for  Avt  resulting  from  station-keeping  analyses  using  Lissajous 
and  halo-type  nominal  orbits.  The  test  for  equivalent  means  between 
populations  of,  say,  halo-type  and  Lissajous  orbit  propellant  costs 

[  39  ] 

(Av^)  requires  the  use  of  a  pooled  variance  equation.  The  pooled 

variance  equation  in  turn  assumes  that  the  sample  variances  are  from 
populations  with  equivalent  variances;  therefore,  the  statistical 
hypothesis  test  for  equal  population  variances  will  logically  preceed 
the  test  for  equal  population  means. 
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A.  Distribution  of  Delta-Velocities 


Histograms  of  Monte  Carlo  simulations  for  the  Delta-Velocity 
Controller  II  and  the  On/Off  Controller  are  shown  in  Figures  3-1  and 
3-2,  respectively.  The  histograms  are  plotted  using  ten  classes  of 
equal  width.  The  vertical  axis  is  the  relative  frequency  of  the  values 
of  Av^  in  that  class.  The  propellant  costs  (values  of  Av^)  for 
station-keeping  on  a  2-year  halo-type  orbit  are  plotted  along  the 
horizontal  axis.  The  class  midpoints  are  labeled;  however,  the 
relative  magnitudes  of  the  Avt  results  for  these  two  controllers  should 
not  be  used  to  compare  the  two  controllers.  This  comparison  of  Av^ 
values  for  the  various  station-keeping  algorithms  will  be  completed 
later  in  this  section.  Both  histograms  show  a  close  resemblance  to  a 
histogram  based  on  the  Gaussian  distribution.  (Recall  that  the  general 
shape  of  a  Gaussian  distribution  is  determined  by  the  variance;  the 
sample  mean  affects  the  expected  frequencies  computed  from  the  Gaussian 
probability  distribution  and  thus  also  helps  to  determine  the  shape  of 
the  histogram  that  could  be  constructed  using  expected  frequencies. ) 
In  Figures  3-1  and  3-2,  the  observed  frequencies  from  the  sample  are 
printed  above  the  classes  in  the  histograms.  The  approximate  expected 
frequencies,  computed  from  the  Gaussian  distribution  using  the  means 
and  standard  deviations  of  the  samples,  are  printed  in  parentheses 
below  the  classes  of  the  histograms. 

The  chi-squared  goodness  of  fit  test  is  appropriately  used  here  to 
verify  mathematically  that  the  frequency  distribution  used  to  construct 
the  given  histogram  fits  the  Gaussian  distribution.  In  order  to 
investigate  this  resemblance,  the  hypotheses  tested  and  the  decision 
rule  are: 


Hypotheses: 

H  distribution  is  Normal, 
o 

H  distribution  is  not  Normal. 

t 

Decision  Rule: 

If  X2  s  ,  conclude  H  . 

o 

Otherwise,  conclude  H  . 

r 
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(.5)  (2.1)  (7.1)  (16)  (25)  (26)  (18)  (7.7)  (2.3)  (.6) 

Total  Delta-Velocities  ( Av  T )  Resulting  from  Controller  II 


Figure  3-1 .  Histogram  for  Total  Delta- Velocities  Resulting  from  use  of  Controller  II. 
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Total  Delta-Velocities  (A  vT )  Resulting  from  On/Off  Controller 


Figure  3-2.  Histogram  of  Total  Delta-Velocities  Resulting  from  use  of  the  On/Off  Controller. 
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Typically,  the  chl-squared  goodness  of  fit  test  includes  the 
listing  of  expected  frequencies  (F^ )  that  are  computed  under  the 

assumption  that  the  null  hypothesis  (Hq)  Is  true.  The  null  hypothesis 
here  Is  that  the  distribution  of  the  values  of  Av  follows  the  Gaussian 

T 

distribution;  therefore,  the  ten  class  boundaries  must  be  determined  in 

terms  of  the  standardized  normal  variable  (Z)  in  preparation  for 

computation  of  the  expected  class  probabilities.  The  class  boundaries 

in  terms  of  the  standardized  normal  variable  can  then  be  used  with  a 

standard  table  of  probability  values  for  the  Gaussian  distribution  to 

compute  expected  probabilities.  The  expected  frequencies  for  each 

class  can  then  be  calculated.  The  expected  frequencies  (F^)  and  the 

observed  frequencies  (the  f  that  can  be  read  directly  from  the 

histograms  or  frequency  distributions)  are  then  used  in  the  computation 

2 

of  the  test  statistic.  The  test  statistic  X  can  be  calculated  as 


k  (f  -  F  ) 

X2  =  £  - - — .  (3-1) 

1=1  F 

l 


Table  3  summarizes  the  mathematics  of  the  chi-squared  goodness  of 
fit  test  for  the  frequency  distribution  of  the  Delta-Velocity 
Controller  II  in  Figure  3-1. 

For  this  test,  the  chi-squared  random  variable  (*2)  has  k-2-1  =  5 
degrees  of  freedom,  reduced  due  to  required  pooling.  (Pooling  is 
required  because  the  expected  frequency  in  any  class  must  be  at  least 
two  for  the  chi-squared  goodness  of  fit  test.  )  The  level  of  confidence 
is  95%.  Therefore,  x  =  11-07.  For  the  distribution  depicted  by  the 
histogram  in  Figure  3-1, 


X 


2 


I 


1  =  1 


«•,-  FI' 


F 

1 


6.639. 
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Table  3.  Chi-Squared  Goodness  of  Fit  Test  for  Figure  3-1. 


Class  Boundaries 

Z  Values 

Probability 

F 

1 

ft 

x2 

l 

.  1415 

.  1681 

-00 

-2.622 

.0044 

•  462  ^ 

\*  . 1233 

.  1645 

.  1946 

-2.622 

-1.973 

.0200 

2. 100 

1 

.1946 

.2211 

-1.973 

-1.325 

.0674 

7.077 

5 

.6096 

.2211 

.2477 

-1 . 325 

-  .676 

.1565 

16.433 

17 

.0196 

.2477 

.2742 

-  .676 

-  .027 

.2397 

25. 169 

31 

1.3512 

.2742 

.3007 

-  .027 

.621 

.2440 

25.620 

20 

1 . 2328 

.3007 

.3273 

.621 

1.278 

.1673 

17.567 

22 

1.1190 

.3273 

.3538 

1.278 

1.919 

.0729 

7.655 

4 

1 . 7447 

.3538 

.3803 

1.919 

2.567 

.0222 

2. 331 . 

3 

[■*  .  4384 

.3803 

.4068 

2.567 

00 

.0052 

.546 

1 

(Note:  The 

8  ymb  o 1 

•  Indicates 

poo  1  1  ng. 

) 

X2=6 . 639 

The  other  frequency  distribution  whose  histogram  is  displayed  in 

2 

Figure  3-2  yields  a  similar  small  value  for  X  .  The  decision  rule  then 

leads  to  the  conclusion  that  the  distribution  of  the  Av  values  is 

T 

normal  for  both  controllers.  For  normally  distributed  random 
variables,  small  sample  sizes  (less  than  30)  can  be  used,  and  there  can 
still  be  confidence  that  the  sample  means  are  normally  distributed. 
Unless  the  distribution  of  the  Avt  values  was  highly  skewed,  a  sample 
size  of  30  would  generally  ensure  that  the  sample  means  were 
approximately  normally  distributed.  Therefore,  confidence  intervals, 
using  the  standard  normal  (Z)  or  student’s  (t)  distribution  can  be 
readily  constructed  using  sample  sizes  of  30. 

One  other  restriction  on  the  minimum  sample  size  could  be  the 
required  width  of  the  confidence  intervals.  Specifically,  the  larger 
the  random  sample  size,  the  narrower  the  resulting  confidence  interval 
for  the  mean.  The  need  for  a  "narrow"  confidence  interval  must 


therefore  be  weighed  against  the  difficulty  in  gathering  the  data.  A 


sample,  as  used  here,  is  a  group  of  station-keeping  (Monte  Carlo) 
simulations  for  which  the  minimal  constraints  are  the  same  and  the 
input  errors  have  identical  means  and  variances.  A  sample  can  then 
also  be  termed  a  "case"  for  later  comparison  with  other  cases 
(samples).  The  sample  size  is  denoted  here  by  the  symbol  "n",  and 
sample  sizes  of  30  are  generally  adequate  to  ensure  that  the  sample 
means  are  at  least  approximately  normally  distributed.  However,  this 
value  of  n  must  also  be  large  enough  to  make  "useful"  confidence 
intervals. 


B.  Sample  Sizes  for  Simulation  Runs 

One  goal  that  motivates  the  use  of  several  random  trials 

(simulations)  of  a  station-keeping  algorithm  is  the  desire  to  construct 

confidence  intervals  for  the  mean  Av  ,  here  denoted  as  A v  for  the 

T 

sample.  A  sample  size  of  perhaps  20  (or  even  40)  for  a  normally 
distributed  random  variable  may  not  provide  a  confidence  interval  that 
is  narrow  enough  for  comparison  purposes.  Too  many  runs  can, 
alternatively,  be  inefficient  while  adding  little  to  the  value  of  the 
study.  The  selection  of  an  appropriate  sample  size  requires  the  choice 
of  confidence  level  (derived  from  the  selected  a  risk)  and  interval 
half  width  (h).  The  standard  deviation  of  the  population  is  the 
parameter  a  that  can  be  estimated  by  the  sample  statistic  s  (sample 
standard  deviation).  The  confidence  interval  for  the  actual  population 
mean  p  is  typically  constructed  surrounding  the  sample  mean  (Av),  such 
that: 


or 


Av  +  (-Z  <r)/ n1/2s  p  s  Av  +  (Z  <r)/n1/2, 


(3-2) 


Av  -h  s  p  s  Av  +  h. 


(3-3) 
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Therefore, 


h  =  (Z  <r)/n1/2. 


(3-4) 


Solving  for  n  results  in  the  following  approximate  guideline  for 
selecting  a  sample  size: 

n  =  (Z  <r)2/h2  a  (Z  s)2/h2.  (3-5) 


Notice  that  this  estimate  for  n  uses  an  approximate  value  (s)  for 
<r  and  a  judgement  concerning  the  size  of  an  acceptable  half  width  (h). 
The  Z  value  is  obtained  from  tables  that  are  available  for  the 
standardized  normal  variable,  using  the  fact  that  one-half  of  the 
a  risk  is  placed  in  each  tail  of  the  distribution.  By  choosing  h  to  be 
less  than  approximately  10%  of  a  Av  value  that  is  equal  to 
.2753,  h  can  be  selected  as  .020  meters  per  second.  The  sample 
standard  deviation  (s)  is  .0409  meters  per  second.  The  standard  normal 
variable,  computed  for  a  =  .01  (a  level  of  confidence  of  99%),  is 
2.576.  Hence, 


n  =  [  (2. 576)  ( .  0409)  ]  2/( .  020) 2  =  27.75  ->  28  (round  up). 


Clearly,  this  is  a  very  approximate  method  to  select  the  sample 
size  n;  however,  sample  sizes  of  30  should  clearly  provide  sufficiently 
accurate  confidence  intervals.  The  value  of  the  half  width  (h)  for  a 
given  confidence  interval  will  depend  on  the  sample’s  standard 
deviation  (s).  In  other  samples,  it  will  undoubtedly  differ  from  the 
.020  meters  per  second  used  in  this  example.  However,  knowing  that  the 
Avt  random  variables  are  approximately  normally  distributed  and,  in  one 
case,  a  very  narrow  confidence  interval  could  be  constructed  with  a 
sample  size  of  28,  further  Monte  Carlo  simulations  are  completed  using 
30  runs  for  each  sample. 
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C.  Results  and  Comparisons 


In  this  section,  results  from  station-keeping  simulations  are 
presented.  Each  station-keeping  method  discussed  earlier  is 

represented.  The  results  are  listed  in  the  order  in  which  the 
algorithms  were  derived,  and  the  population  variances  and  values  of  the 
mean  Avt  (denoted  as  Av)  are  statistically  compared  for  station-keeping 
relative  to  the  halo-type  or  Lissajous  orbit  defined  earlier.  The 
final  subsection  provides  a  survey  of  libration  point  orbit 
station-keeping  costs  listed  in  some  other  works.  Note  that  the 
following  notation  will  be  employed  in  this  section  for  convenience: 
At  is  the  minimum  time  between  control  inputs,  Av  is  the  minimum 

min  min 

control  input  magnitude,  d  is  the  minimal  distance  of  the  spacecraft 

min 

from  the  nominal  path  before  a  control  input  is  allowed  (this  is  then 
also  the  size  of  the  acceptable  torus);  Av  is  an  individual  (vector) 
delta-velocity  input  with  magnitude  Av;  Av^  is  the  (scalar)  total 
delta-velocity  resulting  from  a  single  2-year  station-keeping 
simulation;  "n"  is  the  number  of  simulations  in  a  given  sample  and  is 
often  called  the  sample  size;  Av  is  the  (scalar)  average  of  all  the  Av^ 
values  in  a  given  sample  of  n  simulations;  "s"  is  the  standard 
deviation  from  the  sample  of  n  simulations;  and  a  "case"  is  a  sample  of 
n  station-keeping  simulations  for  which  the  means  and  variances  of  all 
the  input  variables  plus  the  constraints  are  held  at  consistent  levels. 


1.  Delta-Velocity  Controller  I  Results 

The  first  delta  velocity  controller  derived  in  this  work  provides 
excellent  station-keeping  costs  when  the  minimum  time  that  must  elapse 
between  manuevers  (At  )  is  specified  as  40  or  60  days.  The  results 

min 

when  At  is  equal  to  80  days  are  less  promising;  however,  these  costs 

ml  n 

may  be  substantially  reduced  by  an  improved  selection  method  for  the 
target  times  (described  below)  and  weighting  matrix  entries.  Both 
subjects  would  certainly  be  areas  for  future  research  effort. 
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The  target  times  are  chosen  for  this  work  depending  on  the  value 
of  (which  Is  certainly  not  the  only  method  of  choosing  the 

targets)  and  are  summarized  as  one  of  three  possibilities: 

At  =40  days  then  t  =  t  +20  days  and  t  =  t  +40  days, 

At  =60  days  then  t  =  t  +40  days  and  t  =  t  +60  days, 

Bln  10  2  0 

At  =80  days  then  t  =  t  +40  days  and  t  =  t  +80  days. 

“In  1  0  2  0  1 

The  weighting  matrix  entries  are  held  constant  for  these 

comparisons  (other  methods  of  weighting  matrix  computation  could  also 
have  been  used ) : 
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The  results  of  five  samples  (five  cases)  with  sample  sizes  of  30 
each  for  the  2-year  station-keeping  simulations  are  depicted  in 
Table  4.  These  cases  for  Delta-Velocity  Controller  I  are  labeled  in 
the  first  column  as  1-1,  1-2,  etc.  The  column  labeled  as  "At 

min 

indicates  the  minimum  separation  in  days  between  control  inputs  for 

that  case.  The  column  "Av  "  indicates  the  minimum  control  energy  in 

min 

meters  per  second  that  must  be  exceeded  by  a  computed  Av  before  it  will 
be  implemented.  These  minimum  Av  values  agree  with  values  used  in 

(40  41] 

other  libration  point  orbit  station-keeping  studies.  ’  An 

[42] 

investigation  by  Longuski  and  Todd  also  gives  considerable  insight 

into  the  magnitudes  of  the  various  force  levels  affecting  the 

spacecraft.  Their  findings  have  been  used  to  help  determine  the 

minimum  control  energy  levels  that  can  be  used  in  this  effort.  The 

symbol  "d  "  (torus  size)  is  in  kilometers  and  indicates  the  distance 
min 

the  spacecraft  must  be  from  the  nominal  path  before  a  control  force 

will  be  input.  All  three  restrictions  (At  ,  Av  ,  and  d  )  must  be 

min  min  min 

met  or  exceeded  before  a  control  input  is  possible. 

The  average  2-year  Av  for  the  Lissajous  orbit  is  indicated  by 
A  T  J. 

Av  ;  for  the  halo-type  orbit,  it  is  denoted  Av  .  The  sample  standard 
L  h 

deviations  are  s  for  the  Lissajous  orbit  and  s  for  the  halo-type 
orbit.  Sample  means  and  standard  deviations  are  given  in  meters  per 
second.  These  sample  means  and  standard  deviations,  as  listed  in 
Table  4,  are,  in  general,  not  equal  for  the  two  different  types  of 

orbit  within  each  case  (that  is,  Av  *  Av  and  s  *  s  within  each 

L  h  L  H 

case);  however,  the  differences  may  not  be  statistically  significant. 
The  significance  of  the  differences  will  be  discussed  shortly. 

The  results  depicted  in  Table  4  are  derived  from  simulations  using 
the  nominal  orbits  depicted  in  Figures  3-1  and  3-2.  The  nominal  orbits 
are  thought  to  be  representative  of  those  being  considered  for 
near-term  missions.  However,  due  to  the  virtually  infinite  variety  of 
sizes  and  shapes  for  these  types  of  orbits,  no  claim  can  be  made  that 
these  results  would  apply  to  all  libration  point  orbits. 
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Table  4.  Sample  Means  and  Standard  Deviations  for  Controller  I. 


# 

At 

Bln 

(days) 

Av 

aln 

(m/s) 

d 

a  1  n 

(km) 

Lissajous  Orbit 

Halo-Type  Orbit 

Av 

L 

(m/s) 

s 

L 

(m/s) 

Av 

H 

(m/s) 

s 

H 

(m/s) 

1-1 

40 

.015 

0 

.3276 

.0737 

.3395 

.0717 

1-2 

40 

.050 

100 

.7994 

.0903 

.7845 

.  1099 

1-3 

60 

.015 

0 

1.0602 

.7597 

.8945 

.9088 

1-4 

60 

.050 

100 

1 . 3909 

.5180 

1 . 2804 

.  3105 

1-5 

80 

.015 

0 

10.4138 

12.8951 

13.2013 

10.3289 

Table  5  contains  other  useful  data  concerning  the  Delta-Velocity 
Controller  I  samples.  This  data  includes  the  range  of  values  for  the 
random  trials  within  each  case.  For  instance,  in  case  I  —  1 ,  at  least 
one  of  the  random  trials  (one  2-year  Avt)  for  a  Lissajous  orbit  had  the 
minimum  Av^  value  of  .18  meters  per  second.  In  the  same  sample,  at 
least  one  run  had  the  maximum  Av^  value  of  .49  meters  per  second.  The 
range  and  the  standard  deviation  are  both  measures  of  dispersion  for 
the  samples. 


Table  5.  Ranges  for  Av^’s  and  Number  of  Av’ s  for  Controller  I. 


Lissajous  Orbit 

Halo-Type  Orbit 

Average 

# 

At 

nln 

Av 

m  1  n 

Torus 

Av 

L 

Range 

Av 

H 

Range 

Number 

(days) 

(m/s) 

(km) 

(m/s) 

(m/s) 

(m/s) 

(m/s) 

of  Av’s 

1-1 

40 

.015 

0 

.  3276 

. 18-. 49 

.3395 

.23-. 54 

14-15 

1-2 

40 

.050 

100 

.7994 

.65-1.00 

.  7845 

.59-1.03 

10-11 

1-3 

60 

.015 

0 

1 . 0602 

.26-2.92 

.8945 

.27-4.46 

9-10 

1-4 

60 

.050 

100 

1 . 390'.- 

.62-3.02 

1 . 2804 

.74-2.01 

8-9 

1-5 

80 

.015 

0 

10.4138 

.72-31.9 

13.2013 

.63-35.3 

7-9 

The  average  number  of  individual  Av’s  applies  to  both  the 
Lissajous  and  halo-type  orbits;  it  is  the  median  number  of  control 
inputs  required  for  the  2-year  station-keeping  simulations  within  a 
given  case  (sample).  Notice  that  in  case  1-5,  one  or  more  of  the  30 
random  station-keeping  trials  required  only  7  manuevers  over  the  2-year 
period. 

It  is  now  straightforward  to  test  whether  the  "type"  of  orbit 
(Lissajous  or  halo-type)  affects  the  resulting  2-year  station-keeping 
costs.  The  statistical  test  for  equal  means  assumes  approximately 
equal  sample  variances,  hence  that  will  be  a  natural  first  test  to 
conduct.  The  statistical  hypothesis  test  for  equal  variances  assumes 
that  the  test  statistic  (s /s  )  follows  Fisher’s  F  probability 

H  L 

distribution.  The  F  distribution  has  two  types  of  degrees  of  freedom 
(numerator  and  denominator)  and  is  tabulated  also  in  terms  of  the  level 
of  confidence.  The  degrees  of  freedom  are  1  less  than  the  respective 
sample  size  (which  is  n)  for  the  numerator  and  denominator  samples;  so 
that  the  degrees  of  freedom  in  this  example  are  29  (which  is  n-1)  for 
both  the  numerator  and  denominator.  An  F  statistic  denoted  by 
F  indicates  29  degrees  of  freedom  for  both  the  numerator 

(29,29,-995) 

and  denominator  and  a  level  of  confidence  of  99%.  The  a  risk  is  then 

.01  and  is  divided  equally  between  lower  and  upper  bounds  in  the 

*  2  2 

decision  rule.  The  test  statistic  is  F  =  (s/s  ),  and  the  hypotheses 

H  L 

and  decision  rule  are: 


Hypotheses: 

2  2 

HQ:<rH  =  <r^  (Variances  for  halo-type  and  Lissajous  orbits  are  equal.) 

H  :<r2  *  <r2 
1  H  L 


Decision  Rule: 


If  F  s  (s2/s2)  s  F  ,  conclude  H  . 

(29,29,-005)  H  L  (29, 29,. 995)  0 


Otherwise,  conclude  H  . 

l 
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In  the  work  that  follows,  the  lower  F  statistic  in  the  decision  rule  is 
labeled  as  F  ;  the  upper  F  statistic  is  denoted  F  The  F 

low  up 

distribution  value  is  F  _  =  F  =  2.7.  The  value  of  F  is 

(29,29,  .  995)  up  low 

derived  from  the  value  of  F  and  is  given  by  F  =  F  = 

UP  (29,29,-005)  low 

1/2.7  =  .3704. 

The  cases  1-1  through  1-5  in  Table  5  can  now  be  analyzed  for  equal 
variances.  The  numerical  results  and  the  hypothesis  test  conclusions 
are  summarized  in  Table  6. 


Table  6.  Results  of  Equal  Variance  Tests  for  Controller  I. 


# 

s 

L 

(m/s) 

s 

H 

(m/s) 

F 

low 

* 

F 

F  Conclusion 

up 

1-1 

.0737 

.0717 

.3704 

.9465 

2.7 

H 

o 

1-2 

.0903 

.  1099 

.3704 

1.4812 

2.7 

H 

o 

1-3 

.7597 

.9088 

.3704 

1.4310 

2.7 

H 

o 

1-4 

.5180 

.3105 

.3704 

.3593 

2.7 

H 

1 

1-5 

12.8951 

10.3289 

.3704 

.6416 

2.7 

A 

H 

o 


Thus,  in  four  out  of  the  five  cases,  the  variances  can  be  assumed 
to  be  equal.  In  general,  there  is  some  finite  risk  of  finding  in  favor 
of  Hi  even  when  Hq  is  true;  however,  this  a  risk  (.01)  was  chosen  to  be 
very  low  in  this  case.  The  hypothesis  tests  for  equal  population  means 
corresponding  to  case  1-4  should  only  continue  if  some  rough 
equivalence  of  variances  can  be  determined  for  it.  Further  samples  in 
this  case  could  be  compiled  so  that  the  hypothesis  test  for  case  1-4 
could  be  redone,  and  this  could  result  in  substantial  problem  insight. 
However,  the  failure  by  a  small  margin  for  the  hypothesis  test  for  only 
this  case  could  also  allow  the  use  of  a  more  liberal  and  general  rule 
of  thumb:  if  the  sample  variances  are  within  approximately  one  order 
of  magnitude,  the  pooled  variance  equation  can  be  correctly  used. 


74 


Therefore,  the  hypothesis  test  for  equal  population  means  may  be 
pursued  in  all  of  these  cases  1-1  through  1-5. 

The  statistical  hypothesis  test  for  equal  population  means  (pL  for 
Lissajous  and  p  for  halo-type  orbit  delta-velocities)  assumes  that  the 

H 

difference  between  the  two  population  means  is  distributed  according  to 
the  Student’s  t  distribution.  The  test  statistic  is  t  and  is  a 
function  of  the  difference  between  the  sample  means  (Av  -Av  ) . 

H  L 

The  pooled  variance  equation  is 

(n  -1 )s2  +  (n  -1 )s2 

s2=  - b - ■ - E.  (3-40) 

c  (n  -1)  +  (n-1) 

L  H 

where  the  number  of  random  trials  for  both  Lissajous  (nL)  and  halo-type 

(n  )  orbits  is  30  in  this  work.  The  test  statistic  is 
H 


t 


(Av  -Av  ) 

H  L 

s  (1/n  +l/n  )1/2 
c  L  H 


The  hypotheses  and  decision  rule  are: 


Hypotheses: 


H  :p  =  p  . 

0  H 


H  :p  *  p  . 

1  L  H 


Decision  Rule: 

» 

If  -t  s  t  st,  conclude  H  . 

o 

Otherwise,  conclude  H  . 

l 

The  a  risk  of  .01  is  again  equally  divided  in  both  tails  of  the  t 

distribution.  For  this  hypothesis  test,  the  degrees  of  freedom  are 

equal  to  (n  -1)  +  (n  -1),  and  the  t  statistic  is  ±t  =  ±2.66. 

M  L  H  (.995,58) 
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The  numerical  results  and  the  conclusions  are  depicted  in  Table  7. 


Table  7.  Results  of  Equality  of  Means  Tests  for  Controller  I. 


-  A 

A 

• 

# 

Av 

L 

(m/s) 

Av 

H 

(m/s) 

-t 

t 

t 

Conclusion 

1-1 

.3276 

.3395 

-2.66 

.6314 

2.66 

H 

o 

1-2 

.7994 

.7845 

-2.66 

-  5736 

2.66 

H 

o 

1-3 

1 . 0602 

.8945 

-2.66 

-.7661 

2.66 

H 

o 

1-4 

1 . 3909 

1 . 2804 

-2.66 

n/a 

2.66 

n /  a 

1-5 

10.4138 

13.2013 

-2.66 

.9241 

2.66 

H 

o 

The  hypothesis  test  conclusions  are  then,  in  general,  that  there 

is  no  significant  statistical  difference  between  the  2-year 

station-keeping  costs  for  the  Lissajous  and  halo-type  orbits  identified 

in  the  study  when  Delta-Velocity  Controller  I  is  used.  The  entries 

"n/a"  in  Table  7  are  for  the  case  1-4  for  which  the  hypothesis  test 

2  2 

for  equivalent  variances  led  to  the  conclusion  in  favor  of  H  :<r  *  cr  . 

1  H  L 

If  the  variances  were  pooled,  t  for  this  case  would  be  -.9853,  and  the 
conclusion  would  also  be  in  favor  of  equal  means  (H  :u  =  u  ). 

O  L  H 

Now  that  the  statistical  hypothesis  tests  for  equal  population 

variances  and  equal  population  means  between  halo-type  and  Lissajous 

orbits  have  been  completed  for  this  controller,  two  more  short  analyses 

seem  appropriate.  Station-keeping  error  analysis  will  look  at  the 

station-keeping  (propellant)  costs  that  can  be  attributed  to  each  error 

source.  Finally,  confidence  intervals  for  the  mean  Av  cost  for  each 

T 

case  can  be  constructed  here  and  used  later  in  the  comparison 
subsection.  A  limited  station-keeping  error  analysis  for  this 

station-keeping  method  has  been  conducted,  and  the  results  are  shown  in 
Table  8. 

This  error  analysis  seeks  to  quantify  the  relative  contributions 
of  the  individual  error  sources  modeled  in  the  simulations.  The  data 
in  Table  8  corresponds  to  results  for  the  halo-type  orbit  shown  in 
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Figure  1-5.  The  injection  and  tracking  error  levels  are  those  computed 
in  the  chapter  two  orbit  determination  error  analysis  of  this  work  for 
a  halo-type  orbit,  and  the  results  found  here  are  compared  to  those 


Table  8.  Station-Keeping  Error  Analysis  Results. 


( 

Error  Source 

This  Work 

c.  .162) 
Simo 

Injection  Errors 

2.5  % 

n/a 

Control  Errors  (10%) 

2.0% 

10.0% 

Tracking  Errors 

55.9% 

50.0% 

Solar  Reflectivity 

Uncertainty 

35.0% 

15.0% 

Force  Model  and 

Integration 

4.6% 

25.0% 

found  in  Simo.1401  The  major  difference  between  the  error  analysis  in 
this  work  and  the  station-keeping  error  analysis  in  Simo  is  that  his 
study  uses  random  solar  reflectivity  errors  with  a  standard  deviation 
of  2.5%  versus  the  13%  used  here.  A  second  difference  is  that  this 
work  also  includes  gravitational  parameter  uncertainty  in  the  tracking 
error  levels.  (The  nominal  orbit  used  in  Sim6  may  also  be  slightly 
different  from  the  halo-type  orbit  used  here. )  Because  of  the 
increased  level  of  solar  reflectivity  and  gravitational  parameter 
uncertainty,  the  overall  error  level  is  considerably  larger,  and  the 
smaller  contributors  (such  as  control  errors  and  integration  errors) 
provide  a  relatively  smaller  percentage  of  the  total  control  costs. 
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Lastly,  confidence  intervals  for  the  mean  used  in  both  types 
of  orbits  can  also  be  constructed.  The  form  of  these  confidence 
intervals  used  later  in  this  chapter  are 


Lower  Limit  Upper  Limit 


Llssajous 

orbit: 

Av  + 

L 

(-t  SL)/n'L'2* 

s  Av  + 
L 

(t 

\  /  1/2 
s  )/n 

L  L 

Halo-type 

orbit: 

Av  + 

H 

v/nr* "» 

<  Av  + 

H 

(t 

.  1/2 
s  )/n 

H  H 

The  t  statistic  will  have  n  -1  and  n  -1  degrees  of  freedom  for  the 

L  H 

Llssajous  and  halo-type  orbits’  confidence  intervals,  respectively.  An 
a  risk  of  .01  is  used  to  give  a  t  statistic  of  2.462.  The  numerical 
results  are  depicted  in  Table  9. 


Table  9.  Confidence  Intervals  for  Controller  I. 


# 

— iElPBfr  |  BH  !!■!  U  I— 

Lower  Limit 

Upper  Limit 

Lower  Limit 

Upper  Limit 

1-1 

.2946 

.3607 

.3073 

.3717 

1-2 

.7588 

.8400 

.7351 

.  8339 

1-3 

.7187 

1 . 1020 

.4860 

1 . 3030 

1-4 

1.1581 

1 . 6237 

1.1408 

1 . 4200 

1-5 

4.6175 

16.2101 

8.5585 

17.8440 

Notice  that  there  is  considerable  overlap  of  the  intervals  in  each 
case  for  the  halo-type  and  Llssajous  orbits.  These  99%  confidence 
intervals  can  be  used  later  to  compare  the  station-keeping  costs  of 
Delta-Velocity  Controller  1  with  several  other  approaches  such  as  Delta 
Velocity  Controller  II. 
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2.  Delta-Velocity  Controller  II  Results 


The  formulation  of  Delta-Velocity  Controller  II  requires  the 
selection  of  seven  weighting  matrices  (two  more  than  are  needed  for 
controller  I).  The  formulation  also  requires  the  determination  of 
three  target  times.  The  weighting  matrices  used  for  these  simulations 
were  set  for  all  trials  at: 
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The  target  times  are  selected  such  that  if 

At  =60  days  then  t  =t  +  25  days,  t  =t  +  50  days,  and  t  =t  +  75  days, 

or 

At  =80  days  then  t  =t  +  30  days,  t  =t  +  60  days,  and  t  =t  +  90  days. 

This  second  delta-velocity  controller  has  been  evaluated  under 
similar  test  conditions  as  in  cases  1-4  and  1-5  for  the  first 
delta-velocity  controller.  Here,  they  are  labeled  as  cases  1 1 — 1  and 
1 1-2,  respectively.  (See  Table  5  for  parameter  selections  in  cases 

I- 4  and  1-5.  )  Cases  1-4  and  1 1— 1  differ  only  in  the  selection  of 
weighting  matrices  and  target  times;  case  I 1-2  uses  a  larger  Av  and 

"  min 

larger  d  than  case  1-5.  Clearly,  the  same  hypothesis  tests  can  be 

min 

conducted  on  these  controller  II  cases,  and  confidence  intervals  can 
also  be  constructed.  Some  of  the  numerical  results  for  cases  1 1  —  1  and 

II- 2  are  summarized  in  Table  10.  The  notation  used  here  is  identical 
to  that  used  for  Table  4. 


Table  10.  Sample  Means  and  Standard  Deviations  for  Controller  II. 


Lissajous  Orbit 

Halo- 

Type  Orbit 

# 

At 

Av 

d 

Av 

s 

Av 

s 

m  l  n 

■  In 

min 

L 

L 

H 

H 

(days) 

(m/s) 

(km) 

(m/s) 

(m/s) 

(m/s) 

(m/s) 

II-l 

60 

.050 

100 

.8450 

.  1603 

.8124 

.  1233 

1 1-2 

80 

.050 

100 

1 . 0852 

.4682 

1 . 0740 

.  J436 

The  average  Av^  for  case  1-4  was  1.3909  meters  per  second  for  the 
Lissajous  orbit  and  1.2804  mete  per  second  for  the  halo-type  orbit. 
The  improvement  in  station-keeping  costs  is  significant  and  perhaps 
indicates  that  future  work  with  this  controller  may  be  quite  valuable. 
Other  numerical  results  for  Controller  II  are  listed  in  Table  11.  The 
notation  used  here  is  Identical  to  that  used  in  Table  5. 


Table  11.  Sample  Ranges  and  Number  of  Control  Inputs  for  Controller  II. 


Lissajous  Orbit 

Halo-Type  Orbit 

Average 

# 

At 

m  1  n 
(days) 

Av 

min 

(m/s) 

d 

d  1  n 

(km) 

Av 

L 

(m/s) 

Range 

(m/s) 

Av 

H 

(m/s) 

Range 

(m/s) 

Number 

of  Av’ s 

II-l 

60 

.050 

100 

.8450 

.57-1.15 

.8124 

.62-1.08 

7-8 

1 1-2 

80 

.500 

100 

1.0852 

.49-2.40 

1 . 0740 

.63-2.03 

7-8 

Both  cases  show  that  some  simulations  incorporating  this 
station-keeping  strategy  require  as  few  as  seven  manuevers  over  the 
entire  2-year  simulation.  The  range  of  Av^  values  for  both  halo-type 
and  Lissajous  orbits  is  again  listed  as  another  measure  of  the  samples’ 
dispersion.  The  data  in  Tables  10  and  11  can  now  be  tested  for 
equality  of  variance  and  equality  of  population  means  between  halo-type 
and  Lissajous  orbits  within  each  case. 

Initially,  the  hypothesis  that  the  population  variances  are  equal 
is  tested.  The  hypothesis  test  procedures  and  the  notation  are 
Identical  to  those  used  for  these  tests  when  evaluating  results  from 
use  of  Controller  I.  The  hypothesis  test  and  the  decision  rule  are 
summarized  below;  the  numerical  results  and  conclusions  are  listed  in 
Table  12: 


Hypotheses 


H 

o 


<r 


2 

L 


H  : <r 2  *  <r2 

1  H  L 


Decision  Rule 


If 


F 


(29,29,-005) 


s  (s2/s2)  S  F  , 

H  L  (29,29,-995) 


conclude  H  . 

o 


Otherwise,  Conclude  H  . 

l 
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Therefore,  with  99%  confidence,  it  can  be  concluded  that  the 
population  variances  are  equal  when  considering  cases  1 1 — 1  and  1 1-2 
(resulting  from  use  of  the  Delta-Velocity  Controller  II).  Once 
population  variances  are  found  to  be  equal,  the  pooled  variance  can  be 
calculated  and  the  hypothesis  test  for  equal  population  means  can  be 
completed. 


Table  12.  Results  of  Tests  for  Equal  Variances  for  Controller  II. 


« 

s 

L 

(m/s) 

s 

H 

(m/s) 

F 

1  0  w 

* 

F 

F 

up 

Conclusion 

II-l 

.  1603 

.  1233 

.3704 

.5916 

2.7 

H 

o 

1 1-2 

.4682 

.3436 

.3704 

.7339 

2.7 

H 

0 

Again,  the  hypothesis  test  for  equal  population  means  is  conducted 
using  the  t  distribution.  The  hypotheses  and  decision  rule  are 


Hypotheses: 


H  :p  =  p  . 

o  L  H 


H  :p  *  p  . 

1  L  H 


Decision  Rule: 


If  -t  s  t  st,  conclude  H  . 

o 

Otherwise,  conclude  H  . 

l 

The  a  risk  of  .01  Is  again  equally  divided  in  both  tails  of  the  t 

distribution.  The  degrees  of  freedom  are,  for  this  hypothesis  test, 

equal  to  (n  -1)  +  (n  -1)  and  the  t  statistic  is  ±t  =  ±2.66. 

L  H  (.9°5,58) 

The  numerical  results  and  the  conclusions  are  listed  in  Table  13. 
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For  Delta-Velocity  Controller  II,  the  results  clearly  indicate 
that  the  means  for  both  types  of  orbits  for  cases  1 1— 1  and  1 1-2  are 
equal  with  a  99%  level  of  confidence.  Now,  in  order  to  prepare  for 
later  cost  comparisons,  confidence  intervals  for  the  mean  propellant 
costs  can  be  constructed. 

Table  13.  Test  Results  for  Equality  of  Means  for  Controller  II. 


• 

# 

Av 

L 

(m/s) 

Av 

H 

(m/s) 

-t 

t 

t 

Conclusion 

.8450 

.8124 

-2.66 

.4486 

2.66 

H 

0 

1 . 0852 

1 . 0746 

-2.66 

.0344 

2.66 

H 

0 

The  confidence  intervals  within  both  cases  can  be  easily 
formulated,  again  using  the  t  distribution.  An  a  risk  of  .01  is  used 
to  give  a  t  statistic  of  2.462.  The  computed  intervals  for  cases 
that  include  both  Controllers  I  and  II  are  listed  in  Table  14. 


Table 

14.  Confidence 

Intervals  for 

the  Mean  for 

Controller  II. 

# 

Lissajous 

Orbit  (m/s) 

Halo-Type 

Orbit  (m/s) 

Lower  Limit 

Upper  Limit 

Lower  Limit 

Upper  Limit 

1-4 

1.1581 

1 . 6237 

1.1408 

1 . 4200 

II-l 

.7729 

.9171 

.7570 

.8678 

1-5 

4.6175 

16.2101 

8.5585 

17.8440 

1 1-2 

.8747 

1 . 2957 

.9200 

1 . 2285 

Clearly,  Delta-Velocity  Controller  II  exhibits  a  great  improvement 
over  Delta-Velocity  Controller  I  for  cases  involving  At  =  80  days, 

Bln 

that  is,  a  minimum  control  input  separation  time  of  80  days.  The 
Improvement  may  not  completely  be  due  to  the  added  complexity  of 
Controller  II.  In  fact,  Controller  I  may  yield  much  improved  costs 
when  alternative  weighting  matrix  entries  and  target  times  are 
selected.  This  is  certainly  an  area  for  future  inquiry. 
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3.  State  Feedback  Controllers 


The  state  feedback  controller  formulated  in  this  work  can  be  used 

as  a  Continuous  Controller  or  as  an  On/Off  Controller.  The  results 

associated  with  use  of  the  On/Off  Controller  are  primarily  summarized 

here.  The  Continuous  Controller  was  developed  In  early  research  and 

did  not  incorporate  the  realistic  constraint  of  a  minimum  control 

level,  nor  the  use  of  a  minimum  deviation  distance  relative  to  the 

nominal  path.  The  station-keeping  costs  resulting  from  simulations 

that  incorporate  the  Continuous  Controller  are  quite  low;  however, 

spacecraft  thrusters  that  can  deliver  continuous  low  thrust  are  not 

currently  operational,  but  are  in  development  for  use  in  the  next 

decade.  The  low  level  of  thrust  that  may  be  required  to  implement  the 

Continuous  Controller  may  still  be  too  low  for  the  engines  currently  in 

development,  and  the  computed  thrust  level  may  be  of  the  same  order  of 

magnitude  as  the  uncertainty  levels  of  other  forces  on  the  spacecraft. 

As  noted  previously,  in  an  important  work  concerning  this  topic, 

[421 

Longuski  and  Todd  have  done  extensive  research  in  the  area  of 

quantifying  the  various  force  levels  on  a  spacecraft  in  orbit  in  this 
solar  system.  The  results  of  their  work  are  used  to  help  determine  the 
minimal  control  energy  that  Is  practical  for  this  controller. 


a.  Continuous  Controller  Results 

The  costs  resulting  from  application  of  the  Continuous  Controller 
are  briefly  mentioned  here,  and  the  preliminary  results  from  one  Monte 
Carlo  trial  appear  In  Table  15.  These  values  are  computed  for  a 
spacecraft  on  the  2-year  Lissajous  path  (depicted  in  Figure  1-4)  in  the 
vicinity  of  the  interior  libration  point  defined  for  the  Sun-Earth+Moon 
system.  The  error  models  for  each  of  the  six  states  are  listed 
vertically  in  the  same  order  as  the  elements  appear  in  the  residual 
state  vector;  that  is,  x,y,z,x,y,z.  Tracking  and  injection  errors  are 
modeled  at  different  levels  and  are  thus  listed  separately.  The  errors 
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2 

are  represented  in  terms  of  mean  and  variance  in  the  form  N(p,<r  ); 

2 

where,  of  course,  the  mean  is  p  and  the  variance  is  cr  . 

The  error  levels  used  in  Table  15  are  clearly  not  those  listed  in 

Table  1  from  the  orbit  determination  error  analysis  investigation 

[  24 ) 

briefly  described  in  Chapter  1.  This  preliminary  continuous 

control  work  was  completed  well  before  the  error  analysis  was  begun. 
Further  work  on  the  Continuous  Controller  is  not  anticipated  at  this 
time;  the  initial,  partial  results  are  only  presented  here  for 

information  and  comparison  purposes.  The  injection  errors  agree  with 

[32] 

those  found  in  Rodriquez-Canabal  ,  and  the  tracking  errors 

correspond  to  those  of  Sim6[401  for  a  libration  point  orbit. 

Control  input  errors  are  computed  as  a  percent  deviation  from  the 
control  command.  The  control  errors  are  then  input  independently  into 
each  channel.  The  control  Inputs  are  continuous  and  at  a  constant 
magnitude  in  each  control  channel  for  either  30-day  or  45-day  steps  in 
this  study.  The  control  inputs  are  then  recalculated  and  input  at  the 
next  computed  level  for  the  following  30  or  45  days.  This  method  of 
computing  control  Inputs  is  continued  throughout  the  2-year  orbit. 
Position  errors  are  expressed  in  kilometers  and  velocity  errors  in 
meters/second. 

Table  15.  Control  Costs  Associated  with  the  Continuous  Controller. 


Bias  and 

Random  Error 

Sources 

Two-Year  Av 

(meters/second ) 

Injection 

Tracking 

Control 

T 

N( 100,0) 

N(0, 1 . 52} 

N( 100,0) 

N( 100,0) 

N( . 05,0) 

N(0,2.52) 
N(0,152) 
N(0, .0012) 

N  ( 0 , 3% ) 

.  1422 

(30-Day  Steps) 

N( . 05, 0) 

N(0, .0012) 

N( . 05, 0 ) 

N(0, ,0032) 

SAME  AS 

SAME  AS 

SAME  AS 

.5362 

(45-Day  Steps) 

ABOVE 

ABOVE 

ABOVE 

===!^^==== 
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The  constant  acceleration  (In  meters  per  second  per  second)  used 
by  the  controller  over  a  given  period  is  easily  converted  to  a  Av  value 
for  comparison  with  other  controllers.  The  acceleration  is  constant 
over  a  given  time  step,  and  the  equivalent  delta-velocity  can  be 
computed  by  multiplying  the  constant  acceleration  by  the  duration  of 
the  time  step  in  appropriate  units. 

The  method  used  here  provides  results  that  are  approximately  equal 

[33] 

to  those  presented  in  Breakwell  et  al.  In  that  study,  an  integral 

quadratic  cost  function  is  minimized  to  compute  optimal  state  feedback 
control  for  a  periodic  halo-type  orbit  near  the  translunar  libratlon 
point  in  the  Earth-Moon  system.  The  control  accelerations  of  about 
10'8  g’ s  listed  in  Breakwell  et  al  are  slightly  higher  than  the 
8.0xl0  9  g’s  that  are  approximately  required  for  the  method  used  here. 
It  is  interesting  to  note  that  the  solar  reflectivity  uncertainty  at 

-9 

the  two  standard  deviation  level  is  approximately  4.0x10  g’s  for  this 

[42] 

spacecraft  in  the  libratlon  point  orbit.  A  portion  of  the 

difference  in  required  control  levels  may  be  attributed  to  the  much 
larger  value  of  p  in  the  Earth-Moon  system  versus  the  Sun-Earth+Moon 
three-body  problem;  the  orbit  used  in  Breakwell  et  al  is  also  larger, 
in  relation  to  the  relative  size  of  the  respective  three-body  systems, 
than  the  nominal  path  used  here.  In  addition,  Breakwell  et  al  assume 
that  the  only  nonzero  entries  in  the  state  error  weighting  matrix  Q 
(that  appears  in  the  quadratic  cost  function)  are  those  associated  with 
position  errors;  velocity  errors  are  unweighted. 

The  results  listed  in  Table  15  clearly  show  that  the  Continuous 
Controller  can  maintain  the  spacecraft  near  the  nominal  Lissajous  path 
for  2  years  at  a  very  low  cost.  However,  operational  requirements 
suggest  that  the  On/Off  Controller  that  incorporates  the  use  of  some 
minimum  time  between  manuevers,  minimal  control  input  magnitudes,  and  a 
minimum  deviation  distance  from  the  nominal  path  should  be 
investigated.  This  On/Off  Controller,  while  still  using  constant 
accelerative  inputs  over  disjoint  time  periods,  may  more  closely  model 
possible  thrust  devices  being  developed  now  for  use  in  the  next  decade. 
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b.  On/Off  Controller  Results 

The  control  costs  associated  with  the  On/Off  Controller  are  larger 
than  those  found  for  cases  that  incorporate  the  delta-velocity 
controllers  in  this  work.  For  comparison  purposes,  the  accelerative 
inputs  have  been  converted  to  equivalent  delta-velocities  and  have  oeen 
summed  to  reflect  a  2-year  total.  The  Av  value  can  also  be 

Bln 

converted  in  a  similar  way  to  find  minimal  accelerative  force  levels. 
For  instance,  a  minimum  Av  of  .015  meters  per  second  is  equivalent  to  a 

-9 

minimum  constant  acceration  of  approximately  8.681x10  meters  per 

second  per  second  input  for  20  days.  A  minimum  Av  of  .025  meters  per 

second  is  equivalent  to  a  minimum  constant  acceleration  of 

_  8 

approximately  1.4468x10  meters  per  second  per  second.  The  results 
from  simulations  in  two  cases  that  incorporate  the  On/Off  Controller 
are  summarized  in  Table  16.  The  notation  is  identical  to  that  used  in 
previous  sections. 


Table  16.  Mean  and  Variance  Levels  Associated  with  On/Off  Controller. 


Lissajous  Orbit 

Halo-Type 

Orbit 

# 

At 

m  1  n 

(days) 

Av 

m  1  n 

(m/s) 

d 

Bln 

(km) 

Av 

L 

(m/s) 

s 

L 

(m/s) 

Av 

H 

(m/s) 

s 

H 

(m/s) 

III-l 

20 

.015 

100 

.3255 

.0428 

.3661 

.0785 

1 1 1-2 

40 

.025 

100 

.7033 

.2788 

.6999 

.2340 

The 

equivalence  of 

the 

population 

variances 

for  Lissajous  and 

halo-type  orbits  can  now  be  tested.  Again,  Fisher’s  F  test  is  used  to 
test  for  equal  population  variances.  The  a  level  is  set  at  .01,  and 
the  hypotheses  and  decision  rule  are  given  as: 
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Hypotheses : 


H  :  cr 
o 


2 

H 


2 

cr  . 
L 


H  :<yZ  *  a-2. 

l  H  L 
Decision  Rule: 


If 


F 

(29,29,-005) 


S 


(s2/s2)  S  F 

H  L  (29,29,-995) 


conclude  H  . 

0 


Otherwise,  Conclude  H  . 

l 

The  numerical  results  and  conclusions  of  the  hypothesis  tests  for 
equal  variances  are  listed  in  Table  17.  The  notation  is  identical  to 
that  used  in  previous  sections. 

For  both  cases,  the  statistical  hypothesis  test  for  equality  of 

population  variances  within  each  case  for  both  Lissajous  and  halo-type 

orbits  results  in  the  conclusion  that  H  is,  in  fact,  true.  Notice 

o 

that,  if  a  slightly  larger  a  risk  were  used,  cases  1 1 1  —  1  and  1 1 1-2 
would  indicate  that  the  two  populations  do  not  have  equivalent 
variances.  However,  for  these  cases,  we  can,  in  fact,  use  the  pooled 
variance  equation  and  can  then  proceed  with  the  tests  for  equal 
population  means. 


Table  17.  Test  Results  for  Equal  Variances  for  the  On/Off  Controller 


s 

L 

(m/s) 

s 

H 

(m/s) 

F 

I  OH 

* 

F 

F 

up 

Conclusion 

III-l 

.0488 

.0785 

.3704 

2.5876 

2.700 

H 

0 

1 1 1-2 

.2788 

.2340 

.3704 

.7044 

2.700 

H 

0 

The  next  step  in  this  analysis  is  to  test  for  equivalent  Av^ 
population  means.  The  test  again  uses  the  student’s  t  distribution, 
and  an  a  value  of  .01  is  used.  The  hypotheses  and  the  decision  rule 
are  given  by: 
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Hypotheses: 


H  :p  =  n  . 

o  L  H 


H  :n  *  n  . 

l  L 


Decision  Rule: 

* 

If  -t  s  t  st,  conclude  H  . 

o 

Otherwise,  conclude  H  . 

l 


The  a  risk  of  .01  is  again  equally  divided  in  both  tails  of  the  t 

distribution.  The  degrees  of  freedom  are,  for  this  hypothesis  test, 

equal  to  (n  -1)  +  (n  -1)  «»nd  the  t  statistic  is  ±t  =  ±2.66. 

L  H  (.995,58) 

The  numerical  results  and  the  conclusions  of  these  tests  are  listed  in 
Table  18.  The  notation  used  here  is  identical  to  that  used  in  previous 
sections. 


Table  18.  Results  of  Tests  for  Equal  Means  for  the  On/Off  Controller. 


* 

# 

Av 

L 

(m/s) 

Av 

H 

(m/s) 

-t 

t 

t 

Conclusion 

.3255 

.3661 

-2.66 

-2.49 

2.66 

H 

0 

.7033 

.6999 

-2.66 

.1011 

2.66 

H 

0 

For  both  control  options  1 1 1— 1  and  III —2 ,  the  conclusion  is  that 
the  population  means  are  equal  with  99%  confidence.  Clearly,  with  an  a 
risk  slightly  greater  than  1%,  control  option  1 1 1  —  1  would  lead  to  the 
conclusion  that  the  two  populations  have  unequal  means. 

It  might  be  interesting  to  end  this  section  with  a  short 
comparison  of  costs  produced  using  the  On/Off  Controller  with  those 
resulting  from  simulations  with  the  Delta-Velocity  Controller  I.  This 
comparison  should  be  completed  with  the  note  that  the  two  types  of 
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controllers  (delta-velocity  and  on/off)  are,  in  fact,  fundamentally 
different.  Some  cases  are  compared  in  Table  19  by  using  confidence 
intervals  for  the  mean.  The  significant  overlap  of  the  confidence 
Intervals  for  cases  1-1  and  1 1 1 — 1  and  for  cases  1-3  and  1 1 1-2 
illustrate  the  similarity  in  control  costs. 


Table  19.  Comparison  of  Costs  Resulting  from  Use  of  Controller  I  and 
the  On/Off  Controller. 


ft 

Lissajous  Orbit  (m/s) 

Halo  Orbit  (m/s) 

Lower  Limit 

Upper  Limit 

Lower  Limit 

Upper  Limit 

1-1 

.2946 

.3607 

.3073 

.3717 

III-l 

.3063 

.3477 

.3309 

.4013 

1-3 

.7187 

1.1020 

.  4860 

1 . 3030 

1 1 1-2 

.5780 

.8286 

.5947 

.8051 

In  comparing  the  station-keeping  costs,  the  agreement  observed 
between  results  from  use  of  these  two  differing  methods  is  interesting. 
However,  the  On/Off  Controller  produces  costs  that  are  somewhat  higher 
in  general  than  Controller  II,  and  its  performance  when  the  minimum 
manuever  separation  time  is  extended  to  60  or  80  days  is  disappointing. 
Future  work  in  this  area  is  anticipated  and  may  prove  valuable.  The 
last  section  of  this  chapter  contains  a  survey  of  the  station-keeping 
costs  found  in  several  other  libration  point  station-keeping  studies. 
It  also  includes  some  of  the  results  from  this  work. 

4.  Survey  of  Libration  Point  Orbit  Station-Keeping  Costs 

For  completeness,  it  is  important  to  consider  how  libration  point 
orbit  station-keeping  costs  can  be  compared.  Certainly  the  cost  of 
maintaining  the  spacecraft  in  orbit  for  2  years  could  be  a  common 
comparison  value.  However,  the  real  difficulty  here  is  also  closely 
related  to  the  focus  of  the  error  analysis  investigations  in  Chapter  1 
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of  this  work.  A  general  lack  of  consistency  in  input  error  levels  for 
both  orbit  determination  error  analysis  and  station-keeping  studies  for 
libration  point  orbits  has  been  noted.  Some  of  this  inconsistency  can 
be  attributed  to  differing  equipment  and  missions.  Some  cannot.  Two 
papers162’631  concerning  the  same  set  of  studies,  for  instance, 
disagreed  on  the  variances  associated  with  the  tracking  errors  that 
were  used  in  the  station-keeping  simulations. 

Most  studies  do  not  address  the  fact  that  the  control  costs  from  a 
sequence  of  Monte  Carlo  station-keeping  simulations  are  random 
variables.  Each  random  trial  will  provide  a  different  result.  It  is 
assumed  here  that  the  results  published  in  most  studies  are  actually 
the  sample  means  of  the  station-keeping  costs.  One  study  ’  has 
provided  substantial  statistical  details  in  an  excellent 
station-keeping  study  that  could  allow  the  construction  of  confidence 
intervals  for  the  mean  costs  in  that  work.  Such  completeness  of 
data  presentation  is  in  general  lacking  in  other  station-keeping 
studies,  and  any  assumptions  that  must  be  made  about  the  data  will  be 
stated  clearly  in  this  work. 

The  results  of  several  station-keeping  studies  are  presented  in 
Table  20.  The  methods  are  listed  in  order  of  publication  date;  each 
reference  is  listed  by  its  bibliography  number  and  below  it  is  the  date 
of  publication.  Actual  flight  data  is  given  for  ISEE-3.  The  input 
error  levels,  if  available,  are  included  as  standard  deviations  in 
the  order  consistent  with  the  residual  state  vector  (x,y,z;x,y,z) ,  with 
position  state  errors  given  in  kilometers  and  velocity  state  errors 
given  in  milimeters  per  second.  The  standard  deviations  for  solar 
reflectivity  uncertainty  are  listed  as  a  percent  of  the  nominal  value. 
The  control  uncertainty  is  a  one  standard  deviation  error  listed  as  a 
percent  of  the  commanded  control,  usually  input  independently  in  each 
of  the  three  possible  channels. 
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Table  20.  Survey  of  Station-Keeping  Costs  for  Two-Year  Halo  Orbits. 


SUN-EARTH  LI  HALO  ORBIT  STATION -KEEPING  ERRORS 

AND  CONTROL  COSTS 

Solar 

2-Year 

Reflectivity 

Control 

Cost 

Source 

Tracking  Errors  Uncertainty 

Errors 

(m/s) 

Remarks 

(65] 

3,30,30  km; 

N/A 

10% 

15.24 

At  =30  da 

■  In 

(1977) 

15,15,30  mm/s 

26.60 

At  =60  da 

■  In 

I SEE-3 

N/A 

N/A 

N/A 

15.0 

Actual 

[79,80] 

Flight 

(1982) 

Data 

[60] 

2. 7, 3. 9, 3. 4  km; 

10*/. 

N/A 

.08 

4  cm/s/yr 

(1984) 

2 . 4 , 3 . 5 , 1 . 3  mm/s 

[62] 

1.5,2.5,15  km; 

2.5V. 

2.5% 

0.4 

20  cm/s/yr 

(1986) 

1,1,3  mm/s 

[63] 

1.7, 2. 2, 5. 5  km; 

5.0*/. 

2.5% 

0.4 

.7-. 8  m/s/ 

(1987) 

1 . 4, 1 . 4, 2. 4  mm/s 

4  yrs 

[64] 

| (x,y,z)T|s  20  km; 

5.0*/. 

5.0% 

1.3 

4  m/s/ 

(1989) 

| (x,y,z)  |s  15  mm/s 

6  yrs 

[12,91] 

1.5,2.5,15  km; 

N/A 

2.5% 

0.8 

At  =60  da 

■  In 

(1990) 

1,1,3  mm/s 

(one  of 

many  cases) 

This 

1.46,2.64,4.81  km; 

13*/. 

10.0% 

0.81 

At  =60  da 

Work 

1.4,1.85,2.49  mm/s 

( 

.76-. 87) 

a  1  n 

1.1 

At  =80  da 

( 

92-1.23)  "  " 
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The  station-keeping  expenses  are  listed  as  2-year  costs.  When  a 
study  has  printed  the  station-keeping  cost  for  other  than  a  2-year 
value,  this  cost  is  computed  and  the  actual  cost  from  the  study  is 
included  in  the  "remarks"  column.  For  the  control  costs  computed  In 
this  work,  both  the  mean  value  and  then  the  99%  confidence  interval 
limits  (in  parentheses)  are  listed.  If  the  level  of  errors  used  in  a 
given  study  is  not  available,  an  "N/A"  is  indicated  in  the  appropriate 
column.  The  error  levels  used  to  model  orbit  injection  errors  are  not 
listed  in  Table  20.  In  most  of  these  studies,  injection  errors  were 
modeled  at  the  same  level  as  tracking  errors.  In  some  cases,  the 
injection  errors  were  modeled  at  slightly  higher  levels;  however,  at 
these  higher  levels,  the  injection  errors  generally  contribute  only 
about  15%  or  less  to  the  station-keeping  costs. 

Clearly,  there  appears  to  be  a  wide  disparity  in  the 
station-keeping  costs  predicted  for  a  spacecraft  in  a  halo  or  halo-type 
orbit  near  in  the  Sun-Earth  system.  Some  of  these  differences  can 
be  explained  by  the  disagreement  in  error  levels;  some  can  be 
attributed  to  the  variety  of  station-keeping  algorithms.  Also,  these 
studies  do  not  consider  the  same  nominal  orbit.  A  control  cost  such  as 
the  0.8  meters  per  second  for  2  years  found  in  Howell  and 

[9  271 

Pernicka  ’  seems  to  be  quite  reasonable.  This  cost  is  certainly  a 
great  improvement  over  the  7.62  to  13.3  meters  per  second  per  year 
predicted  for  ISEE-3  and  the  actual  mission  results  showing 
approximately  30  meters  per  second  expended  for  less  than  4  years  on 
station.  Further  work  developing  station-keeping  routines  that  can 
tolerate  longer  minimum  times  between  manuevers  (At  ),  increased 

min 

error  levels,  and  the  inclusion  of  a  larger  d  may  be  beneficial. 

min 

Investigations  concerning  the  "optimal"  selection  of  weighting  matrix 
entries  may  enable  existing  methods,  such  as  those  described  in  this 
work,  to  provide  much  improved  results.  Furthermore,  the  selection  of 
target  times  and  control  Input  times  for  these  methods  may  also  be 
fruitful  areas  for  future  research. 
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